Session work scratch: scram X_exit refactor, hot-standby SOS, fat scram tubes, model cheatsheet, journal entry

Multi-session work bundle on a draft branch.  Splits into a clean
sequence of commits later; pushed here so it isn't lost on a reboot.

Reach work
- code/scripts/reach/reach_scram_pj.jl: shutdown_margin halfspace
  X_exit (replaces "n <= 1e-4 AND T_f bound" framing); per-step
  envelope extraction added.
- code/scripts/reach/reach_scram_pj_fat.jl: per-step envelope
  extraction added; shutdown_margin discharge logic mirrored from the
  tight scram script.  3 probes (10/30/60s) all discharge from the
  fat union polytope.
- code/scripts/reach/reach_scram_full_fat.jl (NEW): full nonlinear
  PKE scram reach with fat entry.  Hits the stiffness wall at
  ~1.5 s plant time as expected; saves NaN-tolerant per-step
  envelopes.  Demonstrates concretely why PJ is the right tool for
  the longer-horizon proof.
- code/scripts/reach/reach_heatup_pj.jl: T_REF_START_C constant
  (entry-conditioned ramp) replaces T_STANDBY-init that was making
  the FL controller command cooling at t=0.  Per-step extraction
  already in place.
- code/configs/heatup/tight.toml: bumped maxsteps; probe horizon
  parameterized.

Hot-standby SOS barrier
- code/scripts/barrier/barrier_sos_2d_shutdown.jl (NEW): mirrors the
  operation SOS machinery on the hot-standby thermal projection.
  Includes the eps-slack pattern (so feasibility doesn't silently
  collapse to B == 0).
- code/scripts/barrier/barrier_sos_2d.jl: refactored to use the same
  helper.
- code/src/sos_barrier.jl (NEW): solve_sos_barrier_2d helper module
  factoring out the SOS construction; eps-slack with eps_cap=1.0 to
  avoid unbounded primal.

Library
- code/src/pke_states.jl (NEW): single source of truth for canonical
  initial-condition vectors per DRC mode (op, shutdown, heatup) keyed
  off plant + predicates.
- code/scripts/sim/{main_mode_sweep,validate_pj}.jl, code/CLAUDE.md:
  migrated to pke_states.

Predicates + invariants
- reachability/predicates.json: new shutdown_margin predicate (1%
  dk/k tech-spec floor, expressed as alpha_f*T_f + alpha_c*T_c
  halfspace).  Used as scram X_exit.

Plot script
- code/scripts/plot/plot_reach_tubes.jl: plot_tubes_scram_pj() with
  variant=:fat|:tight knob; plot_tubes_scram_full() for full-PKE
  3-panel (T_c, T_f, rho); plot_tubes_heatup_pj() reads results/
  not reachability/.

Journal + memory
- journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex (NEW):
  long-form entry on the SOS hot-standby barrier and the scram X_exit
  refactor.
- journal/journal.tex: input chain updated.
- claude_memory/ — three new session notes:
  * 2026-04-27-scram-X_exit-shutdown-margin.md
  * 2026-04-28-DICE-2026-conference-intel.md (people, sessions,
    strategic notes for the May 12 talk)
  * 2026-04-28-path1-sos-pj-sketch.md (sketch of nonlinear-SOS via
    polynomial multiply-through; saved for an overnight session)

Docs
- docs/model_cheatsheet.md (NEW): one-page reference of state vector,
  dynamics, constants, modes, predicates, sanity numbers — the talk
  prep cheatsheet Dane asked for.
- docs/figures/reach_*_tubes.png: regenerated with the new mat data.
- presentations/prelim-presentation/outline.md: revised arc per the
  April-28 review pass (cuts: Lyapunov-fails standalone slide,
  operation-tube standalone slide, SOS standalone; adds: scopes-of-
  control framing, scram on the headline result slide).
- app/predicate_explorer.jl: minor.

Hacker-Split: end-of-session scratch bundle

app/predicate_explorer.jl

@@ -1,22 +1,19 @@
 ### A Pluto.jl notebook ###
-# v0.19.40
+# v0.20.24
 
 using Markdown
 using InteractiveUtils
 
-# This Pluto notebook uses @bind for interactivity. The macro is defined locally
+# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
-# so the file remains a valid standalone Julia script when Pluto isn't running.
 macro bind(def, element)
-    #= none:1 =# quote
+    #! format: off
-        local iv = try
+    return quote
-            Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value
+        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
-        catch
-            b -> missing
-        end
         local el = $(esc(element))
         global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
         el
     end
+    #! format: on
 end
 
 # ╔═╡ 9d14e486-1faa-45a9-b235-6367a039f9da
@@ -62,7 +59,9 @@ md"""
 begin
 	pred_path = joinpath(@__DIR__, "..", "reachability", "predicates.json")
     pred_raw  = JSON.parsefile(pred_path)                              
-    md"Loaded `$(relpath(pred_path))`. Top-level keys: $(join(sort(collect(keys(pred_raw))), ", "))."
+    pred_keys_str = join(sort(collect(keys(pred_raw))), ", ")          
+	md"Loaded `$(relpath(pred_path))`. Top-level keys:
+  $(pred_keys_str)."                                                     
 end                                     
 
 # ╔═╡ 48f4bec1-ac48-4318-ba6e-92dbf80a7b2d

claude_memory/2026-04-27-scram-X_exit-shutdown-margin.md

@@ -0,0 +1,88 @@
+# 2026-04-27 — Scram `X_exit` redefinition: shutdown-margin halfspace
+
+## What changed
+
+Replaced the scram-mode `X_exit` predicate from
+`n <= 1e-4 AND T_f <= T_f0 + 50 C` to a single linear halfspace
+`shutdown_margin`:
+
+    alpha_f * T_f + alpha_c * T_c <= -rho_SDM - U_SCRAM + alpha_f*T_f0 + alpha_c*T_c0
+                                  ≈ 0.002972
+
+with `rho_SDM = 0.01` (1% Δk/k tech-spec floor).
+
+Files touched:
+- `reachability/predicates.json` — added `shutdown_margin` under
+  `safety_limits`; updated `mode_definitions.q_scram.X_exit_predicate`
+  and `X_safe_predicate`; left a `_X_exit_history` field for forensics.
+- `code/scripts/reach/reach_scram_pj.jl` — added `RHO_SDM`, `SDM_RHS`
+  constants; reach loop now reports ρ-bounds and the halfspace LHS sup
+  per probe horizon; `.mat` output gets `sdm_lhs_hi`, `rho_max`,
+  `sdm_ok` per horizon plus a global `sdm_rhs`, `rho_sdm`.
+
+## Why
+
+1. **Power threshold was nonlinear in PJ state.** In the prompt-jump
+   reduction, `n = Λ * sum(λ_i*C_i) / (β - ρ)` — ρ depends on `T_f`,
+   `T_c`. So `n ≤ 1e-4` is *not* a halfspace in the reach state. We
+   were reconstructing it post-hoc on every probe, which works for a
+   diagnostic check but doesn't compose with halfspace-based discharge.
+2. **`T_f ≤ T_f0 + 50` was infeasible-by-construction at 60 s.** Decay
+   heat (`Q_sg = 3% P0`) plus a fuel time constant `M_f*C_f / hA ≈ 0.3 s`
+   means `T_f` rapidly equilibrates with `T_c`, but the system loses
+   only ~5 °C in 60 s under constant decay — nowhere near the threshold,
+   and never going back up either. The bound was dressing, not work.
+3. **Shutdown margin is the actual NRC criterion.** Tech specs phrase
+   scram success in Δk/k, not in flux. And ρ is *linear* in
+   `(T_f, T_c)` post-scram (constant `u = U_SCRAM`), so it's a clean
+   single-row halfspace.
+
+## Result
+
+`reach_scram_pj.jl` discharges `shutdown_margin` at all three probe
+horizons (10, 30, 60 s), with massive margin:
+
+| T (s) | reach-sets | wall (s) | ρ at horizon              | discharged |
+|-------|------------|----------|---------------------------|------------|
+| 10    | 6919       | 98.6     | [-0.0507, -0.0504]        | ✓          |
+| 30    | 9900       | 130.5    | [-0.0506, -0.0503]        | ✓          |
+| 60    | 12340      | 164.2    | [-0.0503, -0.0500]        | ✓          |
+
+Required: `ρ ≤ -0.01`. Actual: `|ρ| ≈ 5%`. The Doppler/moderator
+contributions vary by ~3% of `U_SCRAM`, so the margin is dominated by
+rod worth — exactly what physical intuition predicts.
+
+`.mat` output: `results/reach_scram_pj_result.mat`.
+
+## Subtlety I tripped on
+
+Hand-derived `SDM_RHS = 0.00402` using rounded `T_F0=320`, `T_C0=300`.
+Actual values from `pke_params`: `T_C0 = 308.35`, `T_F0 = 328.35`
+(because `DT_CORE = P0/(W_M*C_C) = 36.7 °C`, not 20). The script
+computes `SDM_RHS` from constants so the run was correct (0.002972),
+but the predicate JSON had the stale 0.00402. Fixed by switching the
+`rhs_expr` to a symbolic form. **Rule for next time:** if a constant
+in `predicates.json` is derivable from `pke_params`, write it as a
+symbolic expression, not a baked number — that's what `T_cold0`,
+`T_c0`, `T_standby` are for.
+
+## What's still open
+
+- `X_safe_predicate` still says "fuel_centerline AND cold_leg_subcooled"
+  but the reach script doesn't actually discharge those over the
+  trajectory — only at the probe horizons. Not a problem for the demo
+  (T_f and T_cold are monotone after scram), but the formal obligation
+  is reach-AVOID, not just reach. Worth a follow-up: discharge the
+  invariant halfspaces over the *full flowpipe*, not the endpoint.
+- `prompt_critical_margin_heatup` is the controller-specific PJ-validity
+  predicate. Scram has its own analogous concern (PJ valid only when
+  `β - ρ > 0` with margin). Trivially satisfied here (ρ ≈ -0.05, far
+  from `β = 0.0065`), but worth a parallel `prompt_critical_margin_scram`
+  predicate for completeness — would document the assumption rather than
+  leave it implicit.
+
+## Graduation candidates
+
+- The "rhs_expr should be symbolic when derivable" rule probably belongs
+  in `code/CLAUDE.md` near the predicate-loading section. Hold for now —
+  one occurrence isn't a pattern yet.

claude_memory/2026-04-28-DICE-2026-conference-intel.md

@@ -0,0 +1,268 @@
+# 2026-04-28 — DICE 2026 conference intel (Salt Lake City, May 12-13)
+
+Networking + strategy notes for the 2026 Digital Engineering Conference,
+hosted by INL + University of Utah + Utah Office of Energy Development at
+S.J. Quinney College of Law, U Utah.
+
+## Dane's slot
+
+**Tuesday May 12, 3:30 PM — Breakout Session 10** (afternoon, 2:30–4:30).
+Talk title: *"Leveraging Formal Methods to Build High Assurance Hybrid
+Autonomous Control Systems for Nuclear Power"*. 4th of 6 talks in BS10.
+20-minute slot.
+
+BS10 theme is **risk + assurance**, not tools. Defense-in-depth framing
+(slide 11) lands well here.
+
+## BS10 walkthrough (Dane's session)
+
+| Time | Speaker | Talk | What to know |
+|---|---|---|---|
+| 2:30 | **Olivia Beck** | Metadata Standards for Nuclear Deterrence Test Data | Data-layer / nonproliferation; orthogonal to Dane's work. |
+| 2:50 | **Robert Hayes** | DE of Agility and Risk in Complex SoS | NCSU Associate Professor of Nuclear Engineering. Background is radiation physics, health physics, nonproliferation, dosimetry — *not* his published wheelhouse for SoS/digital engineering. The DICE talk represents a stretch from his usual research; either he's broadening or the topic is a cover for radiation-context work. **Likely Q&A from him: "how does this scale to system-of-systems?"** Have ready: per-mode reach-avoid composition; you verify locally and inherit hybrid correctness. *Possible name collision* — confirm by face on arrival; multiple Robert Hayes in the field. |
+| 3:10 | **Linyu Lin** | Predictive Maintenance Visualization | INL researcher, ML-flavored predictive-maintenance work. Orthogonal but worth a hello — INL collaborator pool. |
+| **3:30** | **Dane** | Formal Methods talk | — |
+| 3:50 | **David Borden** | Cryogenic DT for Neutrino Physics | Specialized; orthogonal. |
+| 4:10 | **Nicole Davis** | "Every Interface is a Risk" | Cyber-flavored closer. Couldn't disambiguate online (very common name). Title strongly suggests OT-cyber posture; she'll be in the same defense-in-depth headspace as Dane's slide 11. **Likely friendly Q&A; mention defense-in-depth out loud and she'll bite.** |
+
+## Top 3 to seek out at the Tuesday reception (4:45–6:30)
+
+### 1. Yue Chen — NREL (likely; common name, see caveat)
+
+**Working hypothesis:** Yue Chen at the National Renewable Energy
+Laboratory (Denver Metro), Ph.D. University of Florida 2012–2016. Coursework
+in optimization/optimal control, stochastic control, control of complex
+networks. Published on combining model-based + model-free methods for
+stochastic control of DERs. NREL has an active **Autonomous Energy
+Systems** thrust with Lyapunov-stability + SOS work for grid-forming
+converters — the talk title (*"Lyapunov-Based Iterative Learning of
+Regions of Attraction for Autonomous Systems"*) fits this lineage.
+
+**Caveat:** "Yue Chen" is extremely common in the field; couldn't 100%
+confirm this is the right Yue Chen. **Verify by face/badge on arrival.**
+LinkedIn: `linkedin.com/in/yuechen10/` for the NREL one.
+
+**Why seek him out:** Methodological neighbor. ROA learning is the
+ML-flavored cousin of Dane's SOS-barrier path. **Conversation opener:**
+*"I'm doing SOS polynomial barriers on a similar problem — what's the
+trade-off in your experience between certified-but-rigid (SOS) and
+learned-but-soft (iterative ROA)?"* Lets him show his work, opens
+collab door.
+
+**If he's at the methodology end (a real stability theorist):** a possible
+collaborator on path 1 (PJ-SOS). Worth a follow-up email after the
+conference.
+
+### 2. Diego Mandelli — Idaho National Laboratory
+
+**Confirmed.** R&D Engineer at INL, Ph.D., works in **Risk Assessment
+and Management Services**. Specializes in dynamic PRA, simulation-based
+risk modeling, AI/data-mining for nuclear safety, knowledge graphs for
+nuclear plant systems. Recent (2024) work on "Technical Language Processing
+of Nuclear Power Plants Equipment Reliability Data" + the MBSE-knowledge-graph
+approach on his DICE talk.
+
+**Why seek him out:** **Licensing/regulatory pathway ally.** Mandelli's
+work is the data + reliability + reasoning layer that complements Dane's
+formal verification — different abstraction levels of the same
+high-assurance problem. He's INL, well-networked, knows the NRC interface.
+
+**Conversation opener:** *"Your KG approach gives a structured reasoning
+layer over reliability data; my work gives bounded-time safety proofs
+over the continuous plant. They're complementary — both feed the
+licensing argument from different sides. How are you seeing the NRC
+respond to formal-methods-based assurance arguments?"*
+
+His DICE talk: BS3 Tuesday morning at 10:25 — *"From Data to Knowledge:
+An MBSE- and Knowledge Graph-Centered AI Framework for Nuclear Reliability
+and Licensing."*
+
+### 3. Sean McBride — Idaho State University
+
+**Confirmed.** Director of the **Informatics Research Institute** at
+Idaho State University's College of Technology. Founded the ICS
+Cybersecurity Associates Degree program at ISU. Background: ex-FireEye
+(built their ICS security business strategy), pioneered DHS ICS-CERT
+threat/vulnerability intelligence, co-founded Critical Intelligence (ICS
+threat intel firm). **One of the people who actually built ICS-CERT.**
+Education-focused now; cares about workforce + students who understand
+both PLCs and physical safeguards.
+
+**Why seek him out:** Closest direct overlap with Dane's formal-methods-
+plus-cyber pitch. McBride represents the OT-cyber audience Dane is trying
+to reach. He's also at INL's neighbor institution — geographic and
+network proximity to Dane's likely collaborators.
+
+His DICE talk: BS6 Tuesday afternoon at 2:30 — *"A Hierarchical Model for
+PLC Code Quality, Safety and Cybersecurity."* Conflicts with Dane's
+slot (BS10), so reception is the moment.
+
+**Conversation opener:** *"Your PLC code quality work and my hybrid
+controller verification work both hit the same target — high-assurance
+control logic — from different layers. I'd love to compare notes on how
+the OT-cyber community is receiving formal-methods arguments."*
+
+## Tuesday morning breakout — strategic room choice
+
+All five rooms (BS1–BS5) run in parallel 9:45–11:45. Pick one and stay.
+Two viable plays:
+
+**Option A — BS1 (methodology overlap):**
+- 10:05 — **Yue Chen**, "Lyapunov-Based Iterative Learning of ROA for Autonomous Systems"
+- 10:45 — Kevin O'Rear (Everstar/Gordian), "DOE→NRC Regulatory Crosswalk via GenAI"
+- 11:05 — Sonali Roy, "AutoDONE: Agentic Framework for NPP Design"
+
+**Option B — BS3 (licensing pathway):**
+- 9:45 — Jieun Lee, "Remote Ops AGN-201 Reactor DT"
+- 10:25 — **Diego Mandelli**, "MBSE + KG AI for Nuclear Reliability and Licensing"
+- 11:05 — **Nicholas Luciano** (ORNL), "Digital Twins to Enable Licensing of Nuclear Innovations"
+
+**Recommendation: Option A.** Yue Chen is the single most valuable
+methodology contact. The licensing-pathway people (Mandelli, Luciano) will
+be at the reception and are findable socially.
+
+## Keynote priorities (12 keynotes; only these 4 matter for Dane)
+
+### Liz Muller — CEO/Co-founder, Deep Fission
+
+**Tuesday 8:45 morning keynote.** Co-founded Deep Fission in 2023 with her
+father Richard Muller (UC Berkeley physicist). Deep Fission deploys
+**off-the-shelf small modular pressurized water reactors a mile underground
+in boreholes** — combining standard PWR tech, deep-borehole drilling
+(oil & gas), and geothermal heat transfer. Reactor named "Gravity."
+Selected for DOE's Reactor Pilot Program; first reactor going in at
+Parsons, Kansas. Recently raised $80M.
+
+Previously co-founded and led Deep Isolation (nuclear waste disposal
+via deep boreholes — same tech family) and Berkeley Earth (climate-data
+nonprofit).
+
+**Why she matters for Dane:** PWRs in unconventional siting → unconventional
+licensing arguments → formal methods becomes more relevant, not less,
+when the regulator can't lean on operational track record. Listen for
+how she frames the regulatory ask. If she emphasizes "we use standard
+PWR tech to minimize licensing risk," that's the opening for
+formal-methods assurance arguments.
+
+### Yasir Arafat — CTO/Co-founder, Aalo Atomics
+
+**Wednesday 9:15 opening keynote.** Founded Westinghouse's eVinci
+microreactor program. Led DOE's **MARVEL project at INL** — first DOE
+reactor authorization in 30 months (very fast). Joined Aalo Atomics from
+INL. Aalo is Austin-based, building the **Aalo-1**: a 10 MWe sodium-cooled
+microreactor inspired by MARVEL, optimized for factory mass-manufacture.
+Partnered with data-center operators. Raised $100M (TechCrunch, Aug 2025).
+
+**Why he matters for Dane:** Formerly at INL (Dane's NRC fellowship
+network), did MARVEL (the real "MARVEL → industry" pipeline). Sodium-
+cooled fast reactor licensing is different from PWR licensing but
+equally hard — autonomous control + formal verification is more, not
+less, valuable. Arafat is the kind of technical founder who'd get
+formal methods immediately.
+
+### Emy Lesofski — Director, Utah Office of Energy Development
+
+Appointed by Gov. Spencer Cox as energy advisor and OED Director in
+late 2023/2024. Previously: U.S. Senate Committee on Appropriations
+(Subcommittee on the Interior, Environment, and Related Agencies) —
+**deep federal appropriations background**. Oversees policy, programs,
+and the Utah San Rafael Energy Lab. UOED has signed an MOU with
+TerraPower exploring siting of an advanced reactor in Utah.
+
+**Why she matters for Dane:** State-level energy authority + federal
+appropriations background = exactly the right node if Dane wants to
+explore state-funded research, advanced-reactor siting work, or the
+San Rafael Energy Lab's research portfolio. UOED is **Diamond sponsor**
+of DICE — she'll be visible and accessible.
+
+### Bryan Lopez — Senior Director, Microsoft Health & Scientific Missions
+
+Federal Strategic Science / Scientific Missions Senior Director, Health
+at Microsoft (Redmond). Previously: DOE Strategic Account Director at
+Microsoft. Earlier: Sandia National Labs, Nuvotech, Air Force Research
+Laboratory. UNM undergrad, U Arizona M.S. (Management Information Systems).
+**He's been the Microsoft↔DOE bridge for years.**
+
+**Why he matters for Dane:** Microsoft is heavy at this conference
+(3 keynote slots — Lopez, Misty Jordan, Nelli Babayan). They have
+discretionary research-engagement budget for federal scientific
+computing. Dane's NRC fellowship + formal methods work is exactly the
+profile Microsoft Federal looks at. Lopez is the contact.
+
+## Other potentially-useful people across the program
+
+- **Nicholas Luciano** (ORNL, BS3 11:05) — R&D in Advanced Engineering
+  Technologies, Nuclear Nonproliferation Division. PhD in nuclear
+  engineering from U Tennessee. Did neutron spectra at SNS, fast-reactor
+  Pu disposition, VVER analysis. Now: digital twins for nuclear
+  licensing — adjacent to Dane.
+- **Max Taylor** (BS2 11:05) — "MBSE to Intrusion Detection Systems."
+  Same defense-in-depth philosophy as Dane.
+- **Prashant Kondle** (BS8 3:50) — "Clearing the Path for AI-Assisted
+  Systems in Regulated Industries." XAI for regulatory acceptance — adjacent
+  to formal methods as a regulatory pathway.
+
+## Hard-question prep (what Dane should expect)
+
+| Person / archetype | Likely question | Prepped answer |
+|---|---|---|
+| **Yue Chen** (or any Lyapunov-ROA person) | "Why SOS over learned ROA? Soundness for adaptivity is a real trade." | Yes, the trade-off is real. Soundness is non-negotiable for NRC. We lose flexibility in exchange for proofs that compose across modes. Complementary, not competitive. |
+| **Robert Hayes** (or any SoS person) | "How does this scale to system-of-systems?" | Per-mode composition. Verify each mode locally; hybrid correctness inherited by composition. Doesn't require monolithic verification. |
+| **Diego Mandelli / SysML-MBSE crowd** | "Why FRET over SysML/MBSE?" | FRET produces machine-checkable LTL; SysML produces human-readable diagrams. Different roles — FRET is downstream of SysML, not a replacement. |
+| **Generative-AI / agentic crowd** (Vaibhav Yadav, Sonali Roy) | "Why not have an LLM do this?" | ML in the safety-critical loop is exactly what we're avoiding. Formal methods give the bounds ML lacks. We're complementary to ML safety analysis, not competitive. **Don't be defensive — the assurance argument is solid.** |
+| **OT-cyber audience** (Sean McBride, Nicole Davis) | "What's your threat model?" | Slide 11 close: formal methods constrain physical-plant behavior even given comms-layer compromise. An assurance axis comms-security alone can't reach. |
+| **Liz Muller / Yasir Arafat archetype** | "How does this help our licensing case?" | Bounded-time safety proofs over the continuous plant give you a quantitative argument for the NRC, not a qualitative one. Verified discrete controller + sound nonlinear reach = "we have proven this can't do the bad thing in this regime." |
+
+## Strategic positioning notes
+
+- **Dane is a methodological outlier.** This conference is heavy on
+  AI/ML/digital-twin/agentic. His formal-methods pitch will stand out —
+  opportunity (memorable) and risk (audience may not be tooled to
+  evaluate it). **Don't apologize.** Lean into the assurance angle;
+  the cyber-leaning subset (BS6, BS7, BS10 second half) gets it
+  instantly.
+- **The regulatory-pathway crowd is the natural ally.** Mandelli,
+  Luciano, O'Rear, Kondle. All asking variants of "how do we get
+  advanced nuclear past the NRC?" Dane has a piece of that puzzle.
+  Find them at the reception.
+- **The Microsoft-Federal triad** (Jordan, Babayan, Lopez) probably
+  has discretionary budget for formal-methods-adjacent federal work.
+  Worth a hello at minimum.
+- **Reception is the highest-leverage window** (Tuesday 4:45–6:30,
+  catered). Wednesday is mostly fireside chats and remarks — less
+  chance to corner the people Dane wants to meet.
+
+## Things to verify on arrival
+
+- **Yue Chen identity** — confirm this is the NREL one (not a different
+  Yue Chen from a different institution). Look at his badge/intro.
+- **Robert Hayes** — confirm whether this is the NCSU radiation-physics
+  professor (research mismatch with talk topic) or a different Robert
+  Hayes. The DICE talk seems out of his published wheelhouse.
+- **Nicole Davis** — couldn't find online; she's almost certainly
+  identifiable from her abstract / intro at the start of her talk
+  in BS10. Decide on the spot whether to ping her after.
+
+## Dane's preferred breakout-session strategy (TL;DR)
+
+1. Tuesday morning 9:45–11:45 → **BS1** (Yue Chen).
+2. Tuesday afternoon 2:30–4:30 → **BS10** (his own session).
+3. Tuesday reception 4:45–6:30 → **find Mandelli + McBride + (Luciano if
+   spotted)**, in that order.
+4. Wednesday morning → keynotes (Arafat); fireside chat is networking-only,
+   doesn't really need close attention.
+5. Anytime he sees Lesofski (Diamond sponsor — she'll be visible) or
+   Lopez — say hi, hand business card.
+
+## Sources
+
+- Diego Mandelli: [INL Researcher Profile](https://bios.inl.gov/Lists/Researcher/DisplayOverrideForm.aspx?ID=538), [Google Scholar](https://scholar.google.com/citations?user=78V6lbsAAAAJ&hl=en)
+- Sean McBride: [ISU Industrial Cybersecurity](https://www.isu.edu/industrialcybersecurity/meet-your-instructors/), [LinkedIn](https://www.linkedin.com/in/sean-mcbride-9705298/)
+- Nicholas Luciano: [ORNL Staff Profile](https://www.ornl.gov/staff-profile/nicholas-p-luciano)
+- Robert Hayes: [NC State Nuclear Engineering](https://ne.ncsu.edu/people/rbhayes/), [Google Scholar](https://scholar.google.com/citations?user=3Jf-ed8AAAAJ&hl=en)
+- Liz Muller: [Deep Fission Leadership](https://www.deepfission.com/about-us/executive-leadership), [LinkedIn](https://www.linkedin.com/in/elizabethmuller/)
+- Yasir Arafat: [Aalo Atomics post](https://www.aalo.com/post/yasir-arafat-of-inls-marvel-to-join-aalo-atomics-as-cto), [LinkedIn](https://www.linkedin.com/in/yasiraalo/)
+- Emy Lesofski: [Cox appointment release](https://governor.utah.gov/press/gov-spencer-cox-appoints-emy-faulkner-lesofski-as-energy-advisor-and-director-of-the-office-of-energy-development/), [LegiStorm bio](https://www.legistorm.com/person/bio/32726/Emelyn_Faulkner_Lesofski.html)
+- Bryan Lopez: [LinkedIn](https://www.linkedin.com/in/bryanlopez/), [ZoomInfo profile](https://www.zoominfo.com/p/Bryan-Lopez/16126421360)
+- Yue Chen (NREL hypothesis): [LinkedIn](https://www.linkedin.com/in/yuechen10/), [NREL Autonomous Energy Systems](https://www.nrel.gov/grid/algorithms)
+- DICE 2026 conference: [INL DICE event page](https://dice.inl.gov/event/digital-engineering-conference-2026/)

claude_memory/2026-04-28-path1-sos-pj-sketch.md

@@ -0,0 +1,125 @@
+# 2026-04-28 — Path 1 sketch: SOS on PJ polynomial dynamics
+
+**Status:** Sketch only. **Do not run yet.** Dane wants this saved for an
+overnight session. Closing the soundness gap on operation/hot-standby
+barriers — moving from linearized SOS to nonlinear-SOS-with-PJ.
+
+## The obstacle
+
+Dynamics in `pke_th_rhs_pj.jl` are *rational*, not polynomial, because
+of the prompt-jump algebraic relation:
+
+    n(x) = Λ · Σ λᵢ Cᵢ / (β - ρ(x))
+
+with `ρ(x) = u + α_f(T_f - T_f0) + α_c(T_c - T_c0)` affine.
+`SumOfSquares.jl` only handles polynomials.
+
+## The trick: multiply through by D(x) = β - ρ(x)
+
+`D(x)` is affine in `(T_f, T_c)` and strictly positive on the operating
+envelope (this is exactly the `prompt_critical_margin_*` predicate —
+already in `predicates.json`).
+
+For each PJ state-equation `dxᵢ/dt = fᵢ(x)/D(x) + polynomial`, multiply:
+
+    D(x) · dxᵢ/dt = fᵢ(x) + D(x) · polynomial
+
+The right-hand side is now polynomial. The Lyapunov/barrier condition
+`dB/dt ≤ 0` becomes (multiplying by `D > 0`):
+
+    D(x) · ∇B(x) · f(x) ≤ 0      ← polynomial inequality
+
+which `SumOfSquares` can handle. Soundness preserved because `D > 0`
+on the operating envelope (and we discharge that as a separate
+predicate, already done in current SOS work via the
+`prompt_critical_margin_*` halfspace).
+
+## Polynomial RHS spelled out
+
+PJ state vector: `x = [C₁..C₆, T_f, T_c, T_cold]`, `u` constant per
+mode. Define `S = Σⱼ λⱼ Cⱼ` (linear in C). Then:
+
+    D · dCᵢ/dt   = βᵢ S - λᵢ Cᵢ D                              ← deg 2
+    D · dT_f/dt  = (P₀ Λ S - hA (T_f - T_c) D) / (M_f c_f)      ← deg 2
+    D · dT_c/dt  = ((hA (T_f - T_c) - 2 W c_c (T_c - T_cold)) D) / (M_c c_c)   ← deg 2
+    D · dT_cold/dt = ((2 W c_c (T_c - T_cold) - Q_sg) D) / (M_sg c_c)           ← deg 2
+
+All right-hand sides are degree ≤ 2 in (C, T_f, T_c, T_cold). With
+B a degree-4 polynomial, the SOS Lie condition `D ∇B · f` is
+degree ≤ 4 + 1 + 2 = 7. Multipliers degree 2. SDP size scales with
+monomial count of the relevant degrees in the relevant variables.
+
+## Dimensionality / tractability
+
+Full 9-state SOS, degree 4: `C(9+4, 4) = 715` monomials. Big SDP
+but tractable on modern CSDP/MOSEK with `chordal sparsity`. May
+need MOSEK rather than CSDP for solver robustness at this scale.
+
+Realistic first attempt: **3D projection on (T_f, T_c, T_cold)** with
+the precursor dynamics treated as bounded uncertainty. Degree 4 in 3
+vars = 35 monomials — same scale as the existing 2D scripts. The
+precursors' contribution to thermal dynamics enters only through
+`S = Σ λⱼ Cⱼ` in the `T_f` equation, and `S` is bounded on the
+reach envelope (we have intervals from prior reach runs). Treat `S`
+as `[S_lo, S_hi]` parametric uncertainty in the SOS — clean Putinar
+multiplier on `[S - S_lo, S_hi - S]`.
+
+This is the pragmatic path: 3D polynomial SOS with bounded `S`.
+Sound for the nonlinear PJ plant on the operating envelope where
+`D > 0` and `S ∈ [S_lo, S_hi]`.
+
+## Session plan for the overnight run
+
+1. Add `D(x)` and `S(x)` symbolic primitives to `sos_barrier.jl`.
+2. Extend `solve_sos_barrier_2d` → `solve_sos_barrier_nd` with
+   parametric-uncertainty multipliers (Putinar on `S ∈ [S_lo, S_hi]`).
+3. Compute `S` envelope from `reach_operation_result.mat` (already
+   exists for operation) and from a long shutdown sim for hot-standby.
+4. Build `barrier_sos_3d_pj_operation.jl`:
+   - 3D state `(δT_f, δT_c, δT_cold)` around operation equilibrium
+   - Polynomial RHS via multiply-through
+   - `D > 0` and `S ∈ [S_lo, S_hi]` as Putinar safety multipliers
+   - Same X_entry, X_unsafe as current 2D operation barrier
+   - Run, compare to linearized result.
+5. Build `barrier_sos_3d_pj_shutdown.jl` — analogous.
+6. If both succeed, write a journal entry making the soundness claim
+   explicit. The new barrier is sound for the PJ-reduced nonlinear
+   plant on `{D > 0} ∩ {S ∈ [S_lo, S_hi]}`. The PJ reduction itself
+   has the Tikhonov error bound (already worked out 2026-04-21).
+   Composing these gives a sound certificate for the *full* nonlinear
+   plant up to the (small, characterized) PJ error.
+
+## Subtleties / things that might bite
+
+- **CSDP may not converge at degree 4 in 3D with multiple Putinar
+  multipliers.** If so, switch to MOSEK (license setup needed) or
+  drop to degree 2 with iterative refinement.
+- **The multiply-through is only valid where D > 0.** If the SOS
+  finds a B that's only valid on a region where the SOS solver's
+  certificate space includes D ≤ 0 points, the result is bogus.
+  Mitigate: use `D` as a Putinar multiplier on the entry/safe sets
+  (so SOS only reasons over `{D > 0}`).
+- **Precursor coupling not certified.** The 3D projection drops the
+  6 precursor states. They're bounded in `S` but their individual
+  dynamics aren't certified. If the SOS proves invariance of
+  `(T_f, T_c, T_cold)` while precursors drift outside `[S_lo, S_hi]`,
+  the certificate is invalid. Need to either (a) certify precursor
+  bounds separately as a forward invariant, or (b) include precursor
+  states in the SOS (back to 9D).
+- **Connection to Tikhonov.** The PJ reduction has error
+  `O(Λ) = O(10⁻⁴)`. Composing PJ-validity-guaranteed barrier + PJ
+  Tikhonov bound = sound certificate for the full plant, modulo the
+  (small) PJ error. Worth working out the full chain rigorously.
+
+## Files that would change
+
+- `code/src/sos_barrier.jl` — generalize to ND + parametric uncertainty
+- `code/scripts/barrier/barrier_sos_3d_pj_operation.jl` — new
+- `code/scripts/barrier/barrier_sos_3d_pj_shutdown.jl` — new
+- `journal/entries/YYYY-MM-DD-pj-sos-soundness.tex` — new
+- `code/CLAUDE.md` — update soundness claim section
+
+## Estimated time
+
+8–12 hours focused work. Realistic overnight run if the SDP is
+well-conditioned. If MOSEK setup eats time, push to a second night.

code/CLAUDE.md

@@ -65,12 +65,15 @@ code/
   README.md                  usage overview
   src/
     pke_params.jl            plant parameters and derived steady state
+    pke_states.jl            canonical ICs for every DRC mode (single source)
     pke_th_rhs.jl            dynamics f(t, x, plant, Q_sg, u)
     pke_linearize.jl         numerical (A, B, B_w) Jacobians
     pke_solver.jl            closed-loop OrdinaryDiffEq driver
     plot_pke_results.jl      4-panel results plot (Plots.jl)
     reach_linear.jl          hand-rolled box reach propagator
     load_predicates.jl       reads reachability/predicates.json
+    sos_barrier.jl           solve_sos_barrier_2d helper (SOS feasibility +
+                             ε-slack pattern; used by barrier_sos_2d{,_shutdown}.jl)
   controllers/
     controllers.jl           ctrl_null, shutdown, heatup, operation,
                              operation_lqr factory, scram
@@ -121,6 +124,15 @@ matches the MATLAB `ode15s` behavior we validated against). Optional
   `Lambda` (~10⁻⁴ s) vs thermal time constants (~10–100 s).
 - **LQR gain is cached** in a `Ref` inside `ctrl_operation_lqr_factory`.
   Reconstruct the factory (not just call it) to retune.
+- **Single-point ICs live in `pke_states(plant; predicates=...)`.** Returns
+  a NamedTuple with `T_standby`, `op`, `shutdown`, `heatup` — every script
+  that needs an IC for forward sim or equilibrium-finding pulls from here.
+  X_entry *polytopes* (for reach) live in `predicates.json` and are loaded
+  separately by reach scripts; don't conflate the two.
+- **SOS barrier programs need ε-slack + scale cap.** The naive feasibility
+  formulation silently returns `B ≡ 0`. `solve_sos_barrier_2d` in
+  `src/sos_barrier.jl` bakes in the (ε > 0, capped) pattern. Use it for
+  any new SOS barrier work.
 - **Plant parameters live in `pke_params()` as a NamedTuple.** For
   reach analysis (TMJets) we duplicate them as `const` globals in
   each reach script — `@taylorize` needs compile-time constants, not

code/configs/heatup/tight.toml

@@ -17,10 +17,21 @@ T_cold_range_C = [278.0, 285.0]
 orderT = 4
 orderQ = 2
 abstol = 1e-9
-maxsteps = 100000
+maxsteps = 1000000
 
 [probes]
-horizons_seconds = [60.0, 300.0]
+# Single probe at the nominal heatup completion time.
+# At T_REF_START_C = 285, RAMP_RATE = 28 C/hr, T_ref reaches T_c0 = 308.35
+# at t = (308.35 - 285) / (28/3600) = 3001 s.  Probing here checks whether
+# the tube has entered X_exit (t_avg_in_range).
+#
+# Why not also probe at T_min = 7714 s as the formal obligation requires:
+# the demo heatup controller has no clamp on T_ref.  Ramping past T_c0
+# would drive the tube past t_avg_high_trip before T_min.  Discharging the
+# full obligation needs a clamped controller (a smooth-min compatible with
+# @taylorize) — flagged as the next thrust.  Tonight we discharge the
+# nominal-heatup-time entry into X_exit; the controller redesign follows.
+horizons_seconds = [3000.0]
 
 [output]
 save_per_step = true

code/scripts/barrier/barrier_sos_2d.jl

@@ -1,27 +1,22 @@
 #!/usr/bin/env julia
 #
-# barrier_sos_2d.jl — SOS polynomial barrier on a 2-state projection.
+# barrier_sos_2d.jl — SOS polynomial barrier on a 2-state projection
+# of the operation-mode LQR closed loop.
 #
 # Proof of concept that SumOfSquares.jl + CSDP can fit a polynomial
-# barrier certificate on a reduced version of the operation-mode
+# barrier certificate on a reduced model.  If this works, scaling to
-# closed-loop.  If this works, scaling to full 10-state is a matter
+# full 10-state is a matter of increasing degree and throughput.
-# of increasing degree and throughput.
 #
-# Reduced dynamics: project the LQR closed-loop onto (dT_c, dn), the
+# Reduced dynamics: project the LQR closed-loop onto (dn, dT_c), the
-# primary safety direction and the dominant unregulated direction.
+# dominant unregulated direction and the primary safety direction.
-# A_red, B_w_red are the 2x2 / 2x1 submatrices corresponding to these
+# A_red, B_red are the 2x2 / 2x1 submatrices corresponding to these
 # components (ignoring cross-coupling into the 8 other states, which is
 # a modeling simplification but keeps the SOS tractable).
 #
 # Safety: |dT_c| ≤ 5 K AND |dn| ≤ 0.15 (i.e. 0.85 ≤ n ≤ 1.15).
 # Entry:  |dT_c| ≤ 0.1 AND |dn| ≤ 0.01.
-# Disturbance: Q_sg deviation |dw| ≤ 0.15·P0.
+# Unsafe focus: dn ≥ +0.15 (high-flux trip; the harder direction
-#
+# under LQR because positive-n excursions trip n_high before T_c trips).
-# Barrier specification (Prajna-Jadbabaie):
-#     B(x) ≤ 0 on X_entry
-#     B(x) ≥ 0 on X_unsafe (= complement of safety)
-#     ∂B/∂x · f(x) ≤ 0 on {B(x) = 0}  (for all w in W)
-# Using SOS multipliers σ_i(x), w-dependence via lossless-disturbance bound.
 
 using Pkg
 Pkg.activate(joinpath(@__DIR__, "..", ".."))
@@ -35,11 +30,11 @@ using CSDP
 include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "sos_barrier.jl"))
 
 plant = pke_params()
-x_op  = pke_initial_conditions(plant)
 
-# Full linearization.
+# Full linearization at full-power steady state.
 A_full, B_full, B_w_full, _, _, _ = pke_linearize(plant)
 
 # Reduced 2x2: rows/cols (1, 9) — n and T_c.
@@ -55,107 +50,46 @@ X_ric, _, _ = arec(A_red, reshape(B_red, :, 1), R_lqr, Matrix(Q_lqr))
 K_red = (R_lqr \ reshape(B_red, 1, :)) * X_ric
 A_cl_red = A_red - reshape(B_red, :, 1) * K_red
 
-println("\n=== SOS barrier attempt — 2-state (n, T_c) projection ===")
+# Cross-coupling check from dropped states.
+cross = A_full[reduce_idx, setdiff(1:10, reduce_idx)]
+
+println("\n=== SOS barrier — 2-state (dn, dT_c) projection of operation LQR ===")
 println("  A_cl_red =")
-show(stdout, "text/plain", A_cl_red)
+show(stdout, "text/plain", A_cl_red); println()
-println()
 println("  B_w_red = $B_w_red")
 println("  eigenvalues: ", round.(eigvals(A_cl_red); sigdigits=4))
+println("  ‖dropped-coupling‖ = $(round(norm(cross); sigdigits=3))")
 println()
 
-# --- SOS formulation ---
+# --- SOS sets ---
-# dx = [dn; dTc] = [x[1]; x[2]] in polynomial variables.
+@polyvar x1 x2   # x1 = dn, x2 = dT_c
-@polyvar x1 x2
 
-# Dynamics with worst-case constant w:
+entry_halfspaces = [
-w_bar = 0.15 * plant.P0
+    0.01 - x1,    # dn ≤ 0.01
-# Split disturbance into its mid + extreme, handle as bounded constant.
+    x1 + 0.01,    # dn ≥ -0.01
-# For the Lie derivative check we use the WORST-CASE w that maximizes
+    0.1  - x2,    # dT_c ≤ 0.1
-# the outward velocity.  Since B_w_red is a known 2-vector and ∂B/∂x
+    x2  + 0.1,    # dT_c ≥ -0.1
-# is polynomial in x, the max-over-w is achieved at w ∈ {-w_bar, +w_bar}.
+]
-# Defer that max — check both worst cases separately.
 
-f_nom = A_cl_red * [x1; x2]            # 2-vector of polynomials in x1, x2
+# Unsafe focus: dn ≥ +0.15 (high-flux trip). Asymmetric — n_high trips
+# at 1.15 (dn = +0.15), n_low at 0.15 (dn = -0.85), so the +0.15
+# direction is the binding one for LQR which has tightly bounded n.
+unsafe_halfspaces = [x1 - 0.15]
 
-# Safety set as intersection of halfspaces g_i ≥ 0:
+# --- Solve ---
-#   g1 = 5 - x2        (dT_c ≤ 5)
+println("  Solving SOS feasibility (degree-4 B, ε-slack capped at 1.0)...")
-#   g2 = x2 + 5        (dT_c ≥ -5)
+result = solve_sos_barrier_2d(A_cl_red, (x1, x2),
-#   g3 = 0.15 - x1     (dn ≤ 0.15)
+                              entry_halfspaces, unsafe_halfspaces;
-#   g4 = x1 + 0.15     (dn ≥ -0.15)
+                              barrier_degree=4, multiplier_degree=2,
-# Unsafe set = complement; for SOS we use the Putinar formulation where
+                              eps_cap=1.0)
-# B ≥ 0 on unsafe.  With multiple unsafe regions (each =complement of
-# one safety halfspace) we'd need one constraint per unsafe region.
-# Simpler: pick one unsafe halfspace to focus on — say n >= 1.15
-# (high-flux trip).  g_u1 = x1 - 0.15.
 
-# Entry set:
+println("  Status: $(result.status)")
-#   g_e1 = 0.1 - x2; g_e2 = x2 + 0.1; g_e3 = 0.01 - x1; g_e4 = x1 + 0.01.
+if result.status == MOI.OPTIMAL && result.ε > 1e-8
-
+    println("  ✅ ε* = $(round(result.ε; digits=4)) — real certificate.")
-g_s1 = 5.0 - x2
+    println("  B(x) = $(result.B)")
-g_s2 = x2 + 5.0
+elseif result.status == MOI.OPTIMAL
-g_s3 = 0.15 - x1
+    println("  ⚠ ε ≈ 0 — solver returned trivial B ≡ 0. No real barrier")
-g_s4 = x1 + 0.15
+    println("     at degree 4 with these sets.")
-g_u_high = x1 - 0.15    # unsafe when n > 1.15 (dn > 0.15)
-g_u_low  = -0.15 - x1   # unsafe when n < 0.85 (dn < -0.15)
-g_e1 = 0.1 - x2
-g_e2 = x2 + 0.1
-g_e3 = 0.01 - x1
-g_e4 = x1 + 0.01
-
-# --- Build the SOS program ---
-solver = optimizer_with_attributes(CSDP.Optimizer, "printlevel" => 0)
-model = SOSModel(solver)
-
-# Barrier polynomial, degree 4.
-monos_B = monomials([x1, x2], 0:4)
-@variable(model, B_poly, Poly(monos_B))
-
-# SOS multipliers for each set constraint, degree 2.
-monos_σ = monomials([x1, x2], 0:2)
-
-# (1) B ≤ 0 on X_entry: -B - Σᵢ σ_eᵢ · g_eᵢ is SOS.
-@variable(model, σ_e1, SOSPoly(monos_σ))
-@variable(model, σ_e2, SOSPoly(monos_σ))
-@variable(model, σ_e3, SOSPoly(monos_σ))
-@variable(model, σ_e4, SOSPoly(monos_σ))
-@constraint(model, -B_poly - σ_e1*g_e1 - σ_e2*g_e2 - σ_e3*g_e3 - σ_e4*g_e4 in SOSCone())
-
-# (2) B ≥ 0 on X_unsafe (using the "high" unsafe region).  Include safety
-# constraints so we stay inside the relevant half:
-#   B - σ_u_high · g_u_high - σ_u_s2 · g_s2 - σ_u_s3 · (-1) is SOS (dummy)
-# Actually: unsafe-high = {x1 ≥ 0.15} alone (unconstrained in x2).
-# Simplest form:
-@variable(model, σ_u, SOSPoly(monos_σ))
-@constraint(model, B_poly - σ_u * g_u_high in SOSCone())
-
-# (3) Lie derivative:  ∇B · f ≤ 0 EVERYWHERE (not just on B=0 boundary).
-# Stronger than needed, but keeps the SDP convex.  The bilinear
-# Putinar form -(∇B·f) - σ_b·B ≥ SOS requires iterative BMI methods;
-# we skip that for this first attempt and use the stronger "global
-# decrease" condition.  If the Hurwitz system admits a quadratic B
-# this should still be solvable.
-dB_dx = [differentiate(B_poly, x1), differentiate(B_poly, x2)]
-# B_w_red is [0, 0] in this projection (Q_sg doesn't directly couple
-# into n or T_c in the linearization), so the disturbance term drops
-# out and the Lie-derivative condition simplifies.
-f_tot = A_cl_red * [x1; x2]
-lie = dB_dx[1] * f_tot[1] + dB_dx[2] * f_tot[2]
-@constraint(model, -lie in SOSCone())
-
-# Feasibility problem — no objective needed.  Any B that satisfies the
-# three SOS constraints is a valid barrier.
-
-println("  Solving SOS program (CSDP)…")
-optimize!(model)
-status = termination_status(model)
-println("  Status: $status")
-if status == MOI.OPTIMAL
-    println("  ✅ SOS barrier found.")
-    println("  B(x) = ", round(value(B_poly); digits=4))
-elseif status == MOI.INFEASIBLE
-    println("  ❌ SOS program infeasible — no degree-4 polynomial B exists")
-    println("     with the given sets and dynamics.  Try higher degree,")
-    println("     larger X_unsafe margin, or different formulation.")
 else
-    println("  ⚠ Solver stopped with: $status")
+    println("  ❌ $(result.status). Try higher degree or relax sets.")
 end

code/scripts/barrier/barrier_sos_2d_shutdown.jl

@@ -0,0 +1,119 @@
+#!/usr/bin/env julia
+#
+# barrier_sos_2d_shutdown.jl — SOS polynomial barrier on the hot-standby
+# (q_shutdown) closed-loop, on a 2-state thermal projection.
+#
+# Companion to barrier_sos_2d.jl (operation-mode LQR projection).
+#
+# Hot-standby controller: u = -5*beta (constant rod insertion). With
+# Q_sg = 0 (no SG load) and small n, the closed loop has a thermal
+# equilibrium where rod-induced negative reactivity balances temperature
+# feedback and decay-heat balances are negligible. We:
+#   1. Find that equilibrium by long-horizon simulation.
+#   2. Linearize there.
+#   3. Reduce to (T_c, T_cold) — the slow safety-relevant thermal modes.
+#      n is decoupled at low power and not the safety driver in this mode.
+#   4. Build a degree-4 SOS barrier on the reduced closed-loop.
+#
+# Safety set (deviation from equilibrium):
+#   |dT_c|    <= 10 K
+#   |dT_cold| <= 15 K
+# Entry set (X_entry from predicates.json::q_shutdown, recentered on x_eq):
+#   |dT_c|    <= 5 K
+#   |dT_cold| <= 5 K
+#
+# X_unsafe-high focus: dT_c >= 10 (over-warming → leaving hot-standby on
+# the wrong side; would imply the shutdown controller cannot hold the
+# plant cool enough).
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
+
+using Printf
+using LinearAlgebra
+using OrdinaryDiffEq
+using DynamicPolynomials
+using SumOfSquares
+using CSDP
+using JSON
+
+include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_states.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_solver.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "sos_barrier.jl"))
+include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
+
+plant = pke_params()
+
+pred_path = joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")
+predicates = JSON.parsefile(pred_path)
+states = pke_states(plant; predicates=predicates)
+println("T_standby = $(round(states.T_standby; digits=2)) °C")
+
+# --- (1) Simulate to quasi-equilibrium ---
+println("\n=== (1) Finding shutdown quasi-equilibrium ===")
+Q_shut = t -> 0.0
+t_sim, X_sim, U_sim = pke_solver(plant, Q_shut, ctrl_shutdown, nothing,
+                                  (0.0, 50000.0); x0=states.shutdown, verbose=false)
+x_eq = X_sim[end, :]
+@printf "  After %.0f s: n=%.3e  T_f=%.2f  T_c=%.2f  T_cold=%.2f\n" t_sim[end] x_eq[1] x_eq[8] x_eq[9] x_eq[10]
+@printf "  Final-step rate ‖dx/dt‖ ≈ %.3e (should be ~0)\n" norm(
+    pke_th_rhs(x_eq, t_sim[end], plant, Q_shut, U_sim[end]))
+
+# --- (2) Linearize at x_eq with u_star = U_SHUTDOWN, Q_star = 0 ---
+println("\n=== (2) Linearizing at x_eq ===")
+u_eq = -5.0 * plant.beta
+A_full, B_full, B_w_full, _, _, _ = pke_linearize(plant;
+    x_star=x_eq, u_star=u_eq, Q_star=0.0)
+
+# u is constant → closed-loop A_cl = A_full (no state feedback term).
+A_cl = A_full
+
+# --- (3) Reduce to (T_c, T_cold) = state indices 9, 10 ---
+reduce_idx = [9, 10]
+A_red = A_cl[reduce_idx, reduce_idx]
+cross = A_cl[reduce_idx, setdiff(1:10, reduce_idx)]
+println("  A_red =")
+show(stdout, "text/plain", A_red); println()
+println("  eigenvalues: ", round.(eigvals(A_red); sigdigits=4))
+println("  ‖dropped-coupling‖ = $(round(norm(cross); sigdigits=3))")
+
+# --- (4) SOS feasibility ---
+println("\n=== (4) SOS barrier (degree-4, ε-slack capped at 1.0) ===")
+@polyvar x1 x2   # x1 = δT_c, x2 = δT_cold
+
+entry_halfspaces = [
+    5.0 - x1,
+    x1 + 5.0,
+    5.0 - x2,
+    x2 + 5.0,
+]
+
+# Unsafe-high focus: dT_c ≥ +10. Asymmetric — over-warming is the harder
+# case for a shutdown controller (rods already maxed in negative; can't
+# add more negative reactivity to compensate).
+unsafe_halfspaces = [x1 - 10.0]
+
+result = solve_sos_barrier_2d(A_red, (x1, x2),
+                              entry_halfspaces, unsafe_halfspaces;
+                              barrier_degree=4, multiplier_degree=2,
+                              eps_cap=1.0)
+
+println("  Status: $(result.status)")
+if result.status == MOI.OPTIMAL && result.ε > 1e-8
+    println(@sprintf "  ✅ ε* = %.4e — real certificate." result.ε)
+    println("  B(x) = $(result.B)")
+elseif result.status == MOI.OPTIMAL
+    println("  ⚠ ε ≈ 0 — solver returned trivial B ≡ 0. No real barrier")
+    println("     at degree 4 for these sets.")
+else
+    println("  ❌ $(result.status). Try higher degree or 3-D extension")
+    println("     (include T_f, since dropped-coupling is non-trivial).")
+end
+
+println("\n=== Equilibrium summary ===")
+@printf "  x_eq = (n=%.3e, T_f=%.3f, T_c=%.3f, T_cold=%.3f)\n" x_eq[1] x_eq[8] x_eq[9] x_eq[10]
+@printf "  u_eq = -5*beta = %.4f\n" u_eq
+@printf "  rho_eq = u_eq + alpha_f*(T_f-T_f0) + alpha_c*(T_c-T_c0) = %.4f\n" (u_eq + plant.alpha_f*(x_eq[8] - plant.T_f0) + plant.alpha_c*(x_eq[9] - plant.T_c0))

code/scripts/plot/plot_reach_tubes.jl

@@ -74,24 +74,35 @@ function plot_tubes_operation()
                  "reach_operation_tubes.png")
 end
 
-function plot_tubes_heatup_pj()
+function plot_tubes_heatup_pj(probe_seconds::Int=8000)
-    mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability",
+    # Reads the per-probe save format from reach_heatup_pj.jl with the
-                       "reach_heatup_pj_tight_full.mat")
+    # latest tight.toml config.  Format: keys prefixed with "T_{probe}_"
+    # for each requested probe horizon (e.g. T_8000_Tc_lo_ts).
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results",
+                       "reach_heatup_pj_tight.mat")
     d = matread(mat_path)
 
-    t_arr   = vec(d["t_arr"])
+    pre = "T_$(probe_seconds)_"
-    Tc_lo   = vec(d["Tc_lo_ts"]);  Tc_hi   = vec(d["Tc_hi_ts"])
+    if !haskey(d, pre * "t_arr")
-    Tf_lo   = vec(d["Tf_lo_ts"]);  Tf_hi   = vec(d["Tf_hi_ts"])
+        # Fall back to legacy bare-key format (older mat files).
-    Tco_lo  = vec(d["Tco_lo_ts"]); Tco_hi  = vec(d["Tco_hi_ts"])
+        pre = ""
-    n_lo    = vec(d["n_lo_ts"]);   n_hi    = vec(d["n_hi_ts"])
+        haskey(d, "t_arr") || error("Neither '$pre' nor bare keys found in $mat_path. Available probes: " *
-    rho_lo  = vec(d["rho_lo_ts"]); rho_hi  = vec(d["rho_hi_ts"])
+            join([replace(k, r"_t_arr$" => "") for k in keys(d) if endswith(k, "_t_arr")], ", "))
+    end
+
+    t_arr   = vec(d[pre * "t_arr"])
+    Tc_lo   = vec(d[pre * "Tc_lo_ts"]);  Tc_hi   = vec(d[pre * "Tc_hi_ts"])
+    Tf_lo   = vec(d[pre * "Tf_lo_ts"]);  Tf_hi   = vec(d[pre * "Tf_hi_ts"])
+    Tco_lo  = vec(d[pre * "Tco_lo_ts"]); Tco_hi  = vec(d[pre * "Tco_hi_ts"])
+    n_lo    = vec(d[pre * "n_lo_ts"]);   n_hi    = vec(d[pre * "n_hi_ts"])
+    rho_lo  = vec(d[pre * "rho_lo_ts"]); rho_hi  = vec(d[pre * "rho_hi_ts"])
 
     Th_lo = 2 .* Tc_lo .- Tco_hi
     Th_hi = 2 .* Tc_hi .- Tco_lo
     dT_lo = 2 .* (Tc_lo .- Tco_hi)
     dT_hi = 2 .* (Tc_hi .- Tco_lo)
 
-    title_stem = "Heatup PJ (tight entry) reach tubes"
+    title_stem = "Heatup PJ (tight entry, T=$probe_seconds s) reach tubes"
     _plot_common(t_arr, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi,
                  dT_lo, dT_hi, rho_lo, rho_hi, n_lo, n_hi, title_stem,
                  "reach_heatup_pj_tubes.png")
@@ -151,18 +162,156 @@ function _plot_common(t, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi,
     println("Saved $outpath")
 end
 
+function plot_tubes_scram_full(probe_seconds::Int=60)
+    # 3-panel scram tube plot using FULL-PKE reach (no prompt-jump
+    # algebraic substitution).  T_c on top, T_f in the middle, ρ on the
+    # bottom.  Per Dane's 2026-04-30 evening request: "I'd love tubes for
+    # t_c t_f and rho."
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results",
+                       "reach_scram_full_fat.mat")
+    d = matread(mat_path)
+    pre = "T_$(probe_seconds)_"
+    haskey(d, pre * "t_arr") || error("No per-step data for probe $probe_seconds in $mat_path")
+
+    t_arr  = vec(d[pre * "t_arr"])
+    Tc_lo  = vec(d[pre * "Tc_lo_ts"]);  Tc_hi  = vec(d[pre * "Tc_hi_ts"])
+    Tf_lo  = vec(d[pre * "Tf_lo_ts"]);  Tf_hi  = vec(d[pre * "Tf_hi_ts"])
+    rho_lo = vec(d[pre * "rho_lo_ts"]); rho_hi = vec(d[pre * "rho_hi_ts"])
+
+    # Convert ρ to dollars.
+    rho_lo_d = rho_lo ./ PLANT.beta
+    rho_hi_d = rho_hi ./ PLANT.beta
+    rho_sdm_dollars = 0.01 / PLANT.beta
+
+    # Panel 1: T_c envelope
+    p1 = plot(xlabel="Time after rod insertion [s]", ylabel="T_c [°C]",
+              title="Average coolant temperature envelope",
+              legend=:right, dpi=180)
+    plot!(p1, t_arr, Tc_hi, fillrange=Tc_lo, fillalpha=0.35,
+          color=:red, linealpha=0, label="T_c tube")
+
+    # Panel 2: T_f envelope
+    p2 = plot(xlabel="Time after rod insertion [s]", ylabel="T_f [°C]",
+              title="Fuel temperature envelope",
+              legend=:right, dpi=180)
+    plot!(p2, t_arr, Tf_hi, fillrange=Tf_lo, fillalpha=0.35,
+          color=:orange, linealpha=0, label="T_f tube")
+
+    # Panel 3: ρ envelope in dollars
+    p3 = plot(xlabel="Time after rod insertion [s]",
+              ylabel="ρ [\$]   (1 \$ = β = prompt-critical)",
+              title="Reactivity envelope",
+              legend=:right, dpi=180)
+    plot!(p3, t_arr, rho_hi_d, fillrange=rho_lo_d, fillalpha=0.35,
+          color=:darkgreen, linealpha=0, label="ρ tube")
+    hline!(p3, [0.0], ls=:dot, color=:black, label="critical")
+    hline!(p3, [-rho_sdm_dollars], ls=:dash, color=:red, lw=2,
+           label="ρ_SDM = -0.01 (\$ = $(round(-rho_sdm_dollars; digits=2)))")
+
+    fig = plot(p1, p2, p3, layout=(3, 1), size=(900, 900),
+               plot_title="Full-PKE scram reach tubes (T=$probe_seconds s, fat entry)")
+    figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
+    isdir(figdir) || mkpath(figdir)
+    outpath = joinpath(figdir, "reach_scram_full_tubes.png")
+    savefig(fig, outpath)
+    println("Saved $outpath")
+end
+
+function plot_tubes_scram_pj(probe_seconds::Int=60; variant::Symbol=:fat)
+    # 2-panel scram tube plot: rho(t) on top, n(t) below.
+    # variant=:tight reads reach_scram_pj_result.mat (small box around
+    # the operating point).  variant=:fat reads reach_scram_pj_fat.mat
+    # (union over shutdown + heatup-tight + operation + LOCA envelopes).
+    mat_filename = variant == :fat ? "reach_scram_pj_fat.mat" :
+                                      "reach_scram_pj_result.mat"
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results", mat_filename)
+    d = matread(mat_path)
+    pre = "T_$(probe_seconds)_"
+    haskey(d, pre * "t_arr") || error("No per-step data for probe $probe_seconds in $mat_path")
+
+    t_arr  = vec(d[pre * "t_arr"])
+    rho_lo = vec(d[pre * "rho_lo_ts"])
+    rho_hi = vec(d[pre * "rho_hi_ts"])
+    n_lo   = vec(d[pre * "n_lo_ts"])
+    n_hi   = vec(d[pre * "n_hi_ts"])
+
+    # Convert rho to dollars (rho / beta).
+    rho_lo_d = rho_lo ./ PLANT.beta
+    rho_hi_d = rho_hi ./ PLANT.beta
+
+    # Panel 1: rho envelope in dollars.
+    p1 = plot(xlabel="Time after rod insertion [s]",
+              ylabel="ρ [\$]   (1 \$ = β = prompt-critical)",
+              title="Reactivity envelope",
+              legend=:right, dpi=180)
+    plot!(p1, t_arr, rho_hi_d, fillrange=rho_lo_d, fillalpha=0.35,
+          color=:darkgreen, linealpha=0, label="ρ tube")
+    hline!(p1, [0.0], ls=:dot, color=:black, label="critical")
+    # Shutdown margin threshold rho_SDM = 0.01 dk/k → in dollars: 0.01/beta
+    rho_sdm_dollars = 0.01 / PLANT.beta
+    hline!(p1, [-rho_sdm_dollars], ls=:dash, color=:red, lw=2,
+           label="ρ_SDM = -0.01 (\$ = $(round(-rho_sdm_dollars; digits=2)))")
+
+    # Panel 2: n envelope (log scale).
+    # Filter NaN / non-positive (PJ reconstruction can produce slack-negatives).
+    n_lo_pos = max.(n_lo, 1e-9)
+    n_hi_pos = max.(n_hi, 1e-9)
+    p2 = plot(xlabel="Time after rod insertion [s]",
+              ylabel="n (normalized power)",
+              title="Power envelope (log scale)",
+              legend=:right, yaxis=:log, dpi=180)
+    plot!(p2, t_arr, n_hi_pos, fillrange=n_lo_pos, fillalpha=0.35,
+          color=:black, linealpha=0, label="n tube")
+    hline!(p2, [1e-4], ls=:dash, color=:orange, lw=2, label="n = 10⁻⁴ (p_above_crit floor)")
+
+    title_suffix = variant == :fat ? "fat entry, T=$probe_seconds s" :
+                                       "T=$probe_seconds s"
+    fig = plot(p1, p2, layout=(2, 1), size=(900, 700),
+               plot_title="Scram reach tubes (PJ-reduced PKE, $title_suffix)")
+    figdir = joinpath(@__DIR__, "..", "..", "..", "docs", "figures")
+    isdir(figdir) || mkpath(figdir)
+    out_filename = variant == :fat ? "reach_scram_pj_fat_tubes.png" :
+                                      "reach_scram_pj_tubes.png"
+    outpath = joinpath(figdir, out_filename)
+    savefig(fig, outpath)
+    println("Saved $outpath")
+end
+
 # CLI dispatch.
 which_plot = length(ARGS) > 0 ? ARGS[1] : "both"
 if which_plot in ("operation", "both")
     plot_tubes_operation()
 end
-if which_plot in ("heatup_pj", "both")
+if which_plot in ("scram_full", "both")
-    mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability",
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results",
-                       "reach_heatup_pj_tight_full.mat")
+                       "reach_scram_full_fat.mat")
+    probe_s = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 60
     if isfile(mat_path)
-        plot_tubes_heatup_pj()
+        plot_tubes_scram_full(probe_s)
     else
-        println("Skipping heatup_pj plot — $mat_path not found.")
+        println("Skipping scram_full plot — $mat_path not found.")
-        println("(Run scripts/reach_heatup_pj_tight_full.jl first.)")
+    end
+end
+if which_plot in ("scram_pj", "both")
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results",
+                       "reach_scram_pj_result.mat")
+    probe_s = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 60
+    if isfile(mat_path)
+        plot_tubes_scram_pj(probe_s)
+    else
+        println("Skipping scram_pj plot — $mat_path not found.")
+        println("(Run scripts/reach/reach_scram_pj.jl first.)")
+    end
+end
+if which_plot in ("heatup_pj", "both")
+    mat_path = joinpath(@__DIR__, "..", "..", "..", "results",
+                       "reach_heatup_pj_tight.mat")
+    # Optional second CLI arg: probe horizon in seconds (default 8000).
+    probe = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 8000
+    if isfile(mat_path)
+        plot_tubes_heatup_pj(probe)
+    else
+        println("Skipping heatup_pj plot — $mat_path not found.")
+        println("(Run scripts/reach/reach_heatup_pj.jl configs/heatup/tight.toml first.)")
     end
 end

code/scripts/reach/reach_heatup_pj.jl

@@ -14,6 +14,12 @@
 #
 # n is algebraic:  n = Λ·Σ λ_i C_i / (β - ρ),  ρ = K_p·(T_ref - T_c).
 #
+# Controller reference: T_ref(t) = T_REF_START_C + RAMP_RATE_CS · t,
+# linear (no clamp at T_c0 — the reach probes when it enters X_exit;
+# clamping would introduce a non-smooth piecewise dynamics @taylorize
+# can't handle).  T_REF_START_C must be at-or-below the X_entry T_c
+# lower bound to keep the controller heating, not cooling.
+#
 # Configuration-driven: pass a TOML config path as the first CLI arg,
 # or omit for the baseline config.
 #
@@ -54,9 +60,18 @@ const T_STANDBY    = T_C0 - 33.333333
 const RAMP_RATE_CS = 28.0 / 3600
 const KP_HEATUP    = 1e-4
 
+# Controller reference starting temperature.  Must be at-or-below the
+# X_entry T_c lower bound so the controller's first move is HEATING, not
+# cooling.  Originally was T_STANDBY = 275, but the heatup X_entry
+# polytope has T_c >= 281 (post-t_avg_above_min), so the FL controller
+# was commanding cooling for the first ~60 s before the ramp caught up.
+# 285 matches `configs/heatup/tight.toml` T_c_lo.  If the entry box
+# changes, update here too.
+const T_REF_START_C = 285.0
+
 # --- Taylorized heatup PJ RHS ---
 @taylorize function rhs_heatup_pj_taylor!(dx, x, p, t)
-    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
+    rho = KP_HEATUP * (T_REF_START_C + RAMP_RATE_CS * x[10] - x[8])
     sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
                 LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
     denom = BETA - rho
@@ -146,8 +161,8 @@ function extract_envelopes(flow_hr)
                    LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
         t_hi_here = high(s, 10)
         t_lo_here = low(s, 10)
-        Tref_lo = T_STANDBY + RAMP_RATE_CS * t_lo_here
+        Tref_lo = T_REF_START_C + RAMP_RATE_CS * t_lo_here
-        Tref_hi = T_STANDBY + RAMP_RATE_CS * t_hi_here
+        Tref_hi = T_REF_START_C + RAMP_RATE_CS * t_hi_here
         rho_lo_here = KP_HEATUP * (Tref_lo - high(s, 8))
         rho_hi_here = KP_HEATUP * (Tref_hi - low(s, 8))
         rho_lo_ts[k] = rho_lo_here

code/scripts/reach/reach_scram_full_fat.jl

@@ -0,0 +1,284 @@
+#!/usr/bin/env julia
+#
+# reach_scram_full_fat.jl — FULL nonlinear PKE scram reach (no PJ).
+#
+# Companion to reach_scram_pj_fat.jl which uses the prompt-jump-reduced
+# 9-state model.  This version keeps n as a state and uses the full
+# 10-state PKE — captures the prompt jump dynamics directly (rather
+# than via the algebraic n substitution).
+#
+# Dane wanted the "real" prompt jump visible in the tubes, even at the
+# cost of stiffness-driven slowness in TMJets.  The penalty: each
+# integration step is bounded by Λ ≈ 10⁻⁴ s.  60 s plant horizon means
+# very many steps and a real risk of over-approximation slack growth.
+#
+# State (11D with augmented time):
+#   x[1]    = n
+#   x[2..7] = C_1..C_6
+#   x[8]    = T_f
+#   x[9]    = T_c
+#   x[10]   = T_cold
+#   x[11]   = t (augmented time)
+#
+# Constant control u = -8β (rod-fall scram).
+# Q_sg = 3% P0 (decay-heat-level sink).
+#
+# Saves per-step envelopes for T_c, T_f, ρ, n at every probe horizon.
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", ".."))
+
+using LinearAlgebra
+using ReachabilityAnalysis, LazySets
+using JSON
+using MAT
+using Printf
+
+# Plant constants — must match pke_params.
+const LAMBDA  = 1e-4
+const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
+    0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
+const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
+const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
+    0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
+
+const P0   = 1e9
+const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
+    50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
+const ALPHA_F, ALPHA_C = -2.5e-5, -1e-4
+
+const T_COLD0 = 290.0
+const DT_CORE = P0 / (W_M * C_C)
+const T_HOT0  = T_COLD0 + DT_CORE
+const T_C0    = (T_HOT0 + T_COLD0) / 2
+const T_F0    = T_C0 + P0 / HA
+
+const U_SCRAM    = -8 * BETA
+const Q_SG_DECAY = 0.03 * P0
+
+const RHO_SDM = 0.01
+const SDM_RHS = -RHO_SDM - U_SCRAM + ALPHA_F*T_F0 + ALPHA_C*T_C0
+
+# Full-PKE Taylorized scram RHS (11 states with augmented t).
+@taylorize function rhs_scram_full_taylor!(dx, x, p, t)
+    rho = U_SCRAM + ALPHA_F * (x[8] - T_F0) + ALPHA_C * (x[9] - T_C0)
+    sum_lam_C = LAM_1*x[2] + LAM_2*x[3] + LAM_3*x[4] +
+                LAM_4*x[5] + LAM_5*x[6] + LAM_6*x[7]
+
+    # n equation — full PKE (the source of stiffness via 1/Λ ≈ 10^4)
+    dx[1] = ((rho - BETA) / LAMBDA) * x[1] + sum_lam_C
+
+    # Precursor balance equations
+    dx[2] = (BETA_1 / LAMBDA) * x[1] - LAM_1 * x[2]
+    dx[3] = (BETA_2 / LAMBDA) * x[1] - LAM_2 * x[3]
+    dx[4] = (BETA_3 / LAMBDA) * x[1] - LAM_3 * x[4]
+    dx[5] = (BETA_4 / LAMBDA) * x[1] - LAM_4 * x[5]
+    dx[6] = (BETA_5 / LAMBDA) * x[1] - LAM_5 * x[6]
+    dx[7] = (BETA_6 / LAMBDA) * x[1] - LAM_6 * x[7]
+
+    # Thermal-hydraulics
+    dx[8]  = (P0 * x[1] - HA * (x[8] - x[9])) / (M_F * C_F)
+    dx[9]  = (HA * (x[8] - x[9]) - 2 * W_M * C_C * (x[9] - x[10])) / (M_C * C_C)
+    dx[10] = (2 * W_M * C_C * (x[9] - x[10]) - Q_SG_DECAY) / (M_SG * C_C)
+
+    dx[11] = one(x[1])
+    return nothing
+end
+
+# --- Build fat X_entry from union of reach envelopes (11D with n) ---
+results_dir = joinpath(@__DIR__, "..", "..", "..", "results")
+pred_raw = JSON.parsefile(joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json"))
+
+# Helper: compute equilibrium C_i for a given n.
+C_eq(n) = [BETA_1/(LAM_1*LAMBDA)*n, BETA_2/(LAM_2*LAMBDA)*n,
+           BETA_3/(LAM_3*LAMBDA)*n, BETA_4/(LAM_4*LAMBDA)*n,
+           BETA_5/(LAM_5*LAMBDA)*n, BETA_6/(LAM_6*LAMBDA)*n]
+
+# 1. Hot-standby (n ≈ 1e-7..1e-4, T near 275 °C).
+sh = pred_raw["mode_boundaries"]["q_shutdown"]["X_entry_polytope"]
+hs_n_lo, hs_n_hi   = sh["n_range"]
+shutdown_lo = [hs_n_lo; C_eq(hs_n_lo); sh["T_f_range_C"][1];
+               sh["T_c_range_C"][1]; sh["T_cold_range_C"][1]]
+shutdown_hi = [hs_n_hi; C_eq(hs_n_hi); sh["T_f_range_C"][2];
+               sh["T_c_range_C"][2]; sh["T_cold_range_C"][2]]
+
+# 2. Operation (n near 1.0, full-power thermals).
+op_path = joinpath(results_dir, "reach_operation_result.mat")
+operation_lo = nothing; operation_hi = nothing
+if isfile(op_path)
+    o = matread(op_path)
+    X_lo_op = o["X_lo"]; X_hi_op = o["X_hi"]   # 10 × M (n at row 1)
+    operation_lo = [minimum(X_lo_op[i, :]) for i in 1:10]
+    operation_hi = [maximum(X_hi_op[i, :]) for i in 1:10]
+end
+
+# X_entry construction note:
+#   We tried widening to a single hyperrectangle covering shutdown
+#   THROUGH high-flux pre-trip (n=1.15, T_c=320, T_f=348).  TMJets
+#   refused — single hyperrectangle that contains both (n=1.15, T_c=320)
+#   AND (n=1e-7, T_c=270) also contains physically-impossible
+#   combinations like (n=1.15, T_c=270).  Over-approximation slack at
+#   t=0 was so wide TMJets failed at the first step.
+#
+#   Right fix: SET SPLITTING — split entry into multiple narrower boxes,
+#   run reach on each separately, take union of tubes.  Multi-hour
+#   build; flagged for follow-up.
+#
+#   For now: union shutdown + "operation transient" box (which itself
+#   spans n=0.85..1.15, the load-follow + high-flux range).  This gives
+#   the widest entry that TMJets can actually solve.  ~7 orders of
+#   magnitude on n; ~50 °C on temperatures.
+
+# Operation-transient box (covers steady + load-follow + pre-trip,
+# physically consistent — high n correlates with high T_f).
+trans_n_lo, trans_n_hi   = 0.85, 1.15
+trans_Tf_lo, trans_Tf_hi = 320.0, 332.0   # ~T_c0 + P/hA at the relevant n
+trans_Tc_lo, trans_Tc_hi = 305.0, 311.0   # tight around T_c0
+trans_Tco_lo, trans_Tco_hi = 285.0, 295.0
+trans_lo = [trans_n_lo; C_eq(trans_n_lo); trans_Tf_lo; trans_Tc_lo; trans_Tco_lo]
+trans_hi = [trans_n_hi; C_eq(trans_n_hi); trans_Tf_hi; trans_Tc_hi; trans_Tco_hi]
+
+# Union: shutdown + transients box.
+fat_lo = copy(shutdown_lo); fat_hi = copy(shutdown_hi)
+for src_pair in [(trans_lo, trans_hi)]
+    lo, hi = src_pair
+    for i in 1:10
+        fat_lo[i] = min(fat_lo[i], lo[i])
+        fat_hi[i] = max(fat_hi[i], hi[i])
+    end
+end
+
+# Defensive clamping (n must be > 0; precursors must be in physical range).
+fat_lo[1] = max(fat_lo[1], 1e-9)
+C_max = [BETA_i/(LAM_i*LAMBDA) * 1.3 for (BETA_i, LAM_i) in
+         zip([BETA_1,BETA_2,BETA_3,BETA_4,BETA_5,BETA_6],
+             [LAM_1,LAM_2,LAM_3,LAM_4,LAM_5,LAM_6])]
+for i in 1:6
+    fat_hi[1+i] = min(fat_hi[1+i], C_max[i])
+    fat_lo[1+i] = max(fat_lo[1+i], 0.0)
+end
+
+# Build the 11D X0 (add augmented time = 0).
+x_lo = vcat(fat_lo, 0.0)
+x_hi = vcat(fat_hi, 0.0)
+X0 = Hyperrectangle(low=x_lo, high=x_hi)
+
+println("\n=== Fat-entry FULL-PKE scram reach (no PJ) ===")
+println("  X_entry (n):       [", round(fat_lo[1]; sigdigits=3), ", ",
+        round(fat_hi[1]; sigdigits=3), "]")
+println("  X_entry (T_f, °C): [", round(fat_lo[8]; digits=2), ", ",
+        round(fat_hi[8]; digits=2), "]")
+println("  X_entry (T_c, °C): [", round(fat_lo[9]; digits=2), ", ",
+        round(fat_hi[9]; digits=2), "]")
+println("  Constant u = -8β = ", round(U_SCRAM; digits=4))
+println("  Probe horizons: 10 s only (full-PKE stiffness wall + slack growth)")
+
+results = Dict{Float64, Any}()
+for T_probe in (10.0,)
+    println("\n--- Probe T = $T_probe s ---")
+    sys  = BlackBoxContinuousSystem(rhs_scram_full_taylor!, 11)
+    prob = InitialValueProblem(sys, X0)
+    try
+        alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=2_000_000)
+        t_start = time()
+        sol = solve(prob; T=T_probe, alg=alg)
+        elapsed = time() - t_start
+        flow = flowpipe(sol)
+        n_sets = length(flow)
+        println("  TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s wall")
+
+        # Per-step overapproximation can produce NaN radii deep into the
+        # run when slack growth blows out (full-PKE scram is stiff and
+        # this is expected after a few seconds of plant time).  Wrap
+        # each step in try/catch so we save the longest contiguous
+        # NaN-free prefix.
+        n_steps_full = length(flow)
+        t_arr_buf     = Float64[]
+        n_lo_buf      = Float64[]; n_hi_buf  = Float64[]
+        rho_lo_buf    = Float64[]; rho_hi_buf = Float64[]
+        Tc_lo_buf     = Float64[]; Tc_hi_buf = Float64[]
+        Tf_lo_buf     = Float64[]; Tf_hi_buf = Float64[]
+        for R in flow
+            try
+                s = set(overapproximate(R, Hyperrectangle))
+                # Skip if any extracted radius is NaN.
+                tval = high(s, 11)
+                nlo  = low(s, 1);  nhi  = high(s, 1)
+                Tflo = low(s, 8);  Tfhi = high(s, 8)
+                Tclo = low(s, 9);  Tchi = high(s, 9)
+                if any(isnan, (tval, nlo, nhi, Tflo, Tfhi, Tclo, Tchi))
+                    break
+                end
+                push!(t_arr_buf,  tval)
+                push!(n_lo_buf,   nlo);  push!(n_hi_buf,   nhi)
+                push!(Tf_lo_buf,  Tflo); push!(Tf_hi_buf,  Tfhi)
+                push!(Tc_lo_buf,  Tclo); push!(Tc_hi_buf,  Tchi)
+                push!(rho_lo_buf, U_SCRAM + ALPHA_F*(Tfhi - T_F0) + ALPHA_C*(Tchi - T_C0))
+                push!(rho_hi_buf, U_SCRAM + ALPHA_F*(Tflo - T_F0) + ALPHA_C*(Tclo - T_C0))
+            catch
+                break
+            end
+        end
+        n_steps = length(t_arr_buf)
+        t_arr     = t_arr_buf
+        n_lo_ts   = n_lo_buf;   n_hi_ts   = n_hi_buf
+        rho_lo_ts = rho_lo_buf; rho_hi_ts = rho_hi_buf
+        Tc_lo_ts  = Tc_lo_buf;  Tc_hi_ts  = Tc_hi_buf
+        Tf_lo_ts  = Tf_lo_buf;  Tf_hi_ts  = Tf_hi_buf
+
+        # Build a "last" hyperrectangle from the final NaN-free reach set.
+        if n_steps == 0
+            error("All reach-sets produced NaN under overapproximate.")
+        end
+        # Use last buffered values for the "at T_probe" report.
+        last_n_lo = n_lo_ts[end];  last_n_hi = n_hi_ts[end]
+        last_Tf_lo = Tf_lo_ts[end]; last_Tf_hi = Tf_hi_ts[end]
+        last_Tc_lo = Tc_lo_ts[end]; last_Tc_hi = Tc_hi_ts[end]
+        last_rho_lo = rho_lo_ts[end]; last_rho_hi = rho_hi_ts[end]
+        println("  Extracted $n_steps NaN-free reach-sets out of $(n_steps_full) total")
+
+        @printf "  n at last good step:   [%.4e, %.4e]\n" last_n_lo last_n_hi
+        @printf "  T_f at last good step: [%.2f, %.2f] °C\n" last_Tf_lo last_Tf_hi
+        @printf "  T_c at last good step: [%.2f, %.2f] °C\n" last_Tc_lo last_Tc_hi
+        @printf "  ρ at last good step:   [%.5f, %.5f]\n" last_rho_lo last_rho_hi
+        @printf "  last good time: %.4f s\n" t_arr[end]
+
+        results[T_probe] = (status="OK", n_sets=n_sets, elapsed=elapsed,
+                            t_arr=t_arr,
+                            n_lo_ts=n_lo_ts, n_hi_ts=n_hi_ts,
+                            rho_lo_ts=rho_lo_ts, rho_hi_ts=rho_hi_ts,
+                            Tc_lo_ts=Tc_lo_ts, Tc_hi_ts=Tc_hi_ts,
+                            Tf_lo_ts=Tf_lo_ts, Tf_hi_ts=Tf_hi_ts)
+    catch err
+        msg = sprint(showerror, err)
+        println("  FAILED: ", first(msg, 300))
+        results[T_probe] = (status="FAILED", err=first(msg, 300))
+        break
+    end
+end
+
+# Save .mat — same key naming convention as the PJ scram script.
+mat_out = joinpath(results_dir, "reach_scram_full_fat.mat")
+saved = Dict{String, Any}()
+saved["fat_lo"]   = fat_lo
+saved["fat_hi"]   = fat_hi
+saved["sources"]  = ["shutdown", "operation", "loca"]
+saved["sdm_rhs"]  = SDM_RHS
+saved["rho_sdm"]  = RHO_SDM
+for T_probe in (10.0, 30.0, 60.0)
+    haskey(results, T_probe) || continue
+    r = results[T_probe]
+    r.status == "OK" || continue
+    pre = "T_$(Int(T_probe))_"
+    saved[pre * "t_arr"]     = r.t_arr
+    saved[pre * "n_lo_ts"]   = r.n_lo_ts
+    saved[pre * "n_hi_ts"]   = r.n_hi_ts
+    saved[pre * "rho_lo_ts"] = r.rho_lo_ts
+    saved[pre * "rho_hi_ts"] = r.rho_hi_ts
+    saved[pre * "Tc_lo_ts"]  = r.Tc_lo_ts
+    saved[pre * "Tc_hi_ts"]  = r.Tc_hi_ts
+    saved[pre * "Tf_lo_ts"]  = r.Tf_lo_ts
+    saved[pre * "Tf_hi_ts"]  = r.Tf_hi_ts
+end
+matwrite(mat_out, saved)
+println("\nSaved fat-entry full-PKE scram reach to $mat_out")

code/scripts/reach/reach_scram_pj.jl

@@ -2,9 +2,15 @@
 #
 # reach_scram_pj.jl — nonlinear reach on scram, prompt-jump model.
 #
-# Scram obligation: from any operating-envelope state, drive n down to
+# Scram obligation: from any operating-envelope state, drive total
-# n <= 1e-4 within T_max = 60 s.  Constant control u = -8*beta (rods
+# reactivity below the shutdown-margin threshold (rho <= -0.01, i.e.
-# slammed in).  Q_sg = 3% P0 (decay-heat-level sink, placeholder).
+# 1% dk/k subcritical) within T_max = 60 s.  Constant control
+# u = -8*beta (rods slammed in).  Q_sg = 3% P0 (decay-heat-level sink,
+# placeholder).
+#
+# X_exit halfspace (from reachability/predicates.json::shutdown_margin):
+#   alpha_f * T_f + alpha_c * T_c <= 0.00402
+# discharged when sup over reach set of LHS <= 0.00402.
 #
 # 9-state PJ model (10D with augmented time).
 
@@ -38,6 +44,10 @@ const T_F0    = T_C0 + P0 / HA
 const U_SCRAM = -8 * BETA          # rod worth applied at scram
 const Q_SG_DECAY = 0.03 * P0       # constant decay-heat-level sink
 
+# X_exit threshold: shutdown_margin halfspace, mirrors predicates.json.
+const RHO_SDM = 0.01                                            # 1% dk/k
+const SDM_RHS = -RHO_SDM - U_SCRAM + ALPHA_F*T_F0 + ALPHA_C*T_C0   # ≈ 0.00402
+
 # Taylorized scram RHS, PJ form.
 @taylorize function rhs_scram_pj_taylor!(dx, x, p, t)
     rho = U_SCRAM + ALPHA_F * (x[7] - T_F0) + ALPHA_C * (x[8] - T_C0)
@@ -82,6 +92,7 @@ println("  X_entry: small box around operating point (n ≈ 1.0)")
 println("  Constant u = -8*beta = $(round(U_SCRAM; digits=4))")
 println("  Q_sg = 3% P0 (decay-heat sink)")
 println("  T_max = 60 s")
+println("  X_exit: alpha_f*T_f + alpha_c*T_c <= $(round(SDM_RHS; sigdigits=4))  (rho <= -$(RHO_SDM))")
 
 results = Dict{Float64, Any}()
 for T_probe in (10.0, 30.0, 60.0)
@@ -98,6 +109,33 @@ for T_probe in (10.0, 30.0, 60.0)
         println("  TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s wall")
 
         flow_hr = overapproximate(flow, Hyperrectangle)
+
+        # --- Per-step envelopes for plotting tubes ---
+        n_steps = length(flow_hr)
+        t_arr     = zeros(n_steps)
+        n_lo_ts   = zeros(n_steps); n_hi_ts   = zeros(n_steps)
+        rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps)
+        Tc_lo_ts  = zeros(n_steps); Tc_hi_ts  = zeros(n_steps)
+        Tf_lo_ts  = zeros(n_steps); Tf_hi_ts  = zeros(n_steps)
+        for (k, R) in enumerate(flow_hr)
+            s = set(R)
+            t_arr[k] = high(s, 10)
+            sumLC_lo_k = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
+                         LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
+            sumLC_hi_k = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
+                         LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
+            rho_lo_k = U_SCRAM + ALPHA_F*(high(s,7) - T_F0) + ALPHA_C*(high(s,8) - T_C0)
+            rho_hi_k = U_SCRAM + ALPHA_F*(low(s,7)  - T_F0) + ALPHA_C*(low(s,8)  - T_C0)
+            denom_lo_k = BETA - rho_hi_k
+            denom_hi_k = BETA - rho_lo_k
+            n_lo_ts[k]   = denom_lo_k > 0 ? LAMBDA * sumLC_lo_k / denom_hi_k : NaN
+            n_hi_ts[k]   = denom_lo_k > 0 ? LAMBDA * sumLC_hi_k / denom_lo_k : NaN
+            rho_lo_ts[k] = rho_lo_k
+            rho_hi_ts[k] = rho_hi_k
+            Tc_lo_ts[k]  = low(s, 8);  Tc_hi_ts[k] = high(s, 8)
+            Tf_lo_ts[k]  = low(s, 7);  Tf_hi_ts[k] = high(s, 7)
+        end
+
         # Reconstruct n at last time step from C and T_c.
         last = set(flow_hr[end])
         sumLC_lo = LAM_1*low(last,1) + LAM_2*low(last,2) + LAM_3*low(last,3) +
@@ -113,13 +151,33 @@ for T_probe in (10.0, 30.0, 60.0)
         Tc_final = (low(last, 8), high(last, 8))
         Tf_final = (low(last, 7), high(last, 7))
         Tcold_final = (low(last, 9), high(last, 9))
+        # shutdown_margin halfspace LHS: alpha_f*T_f + alpha_c*T_c.
+        # Coefficients negative → sup over the box at low(T_f), low(T_c).
+        sdm_lhs_hi = ALPHA_F*low(last,7) + ALPHA_C*low(last,8)
+        sdm_lhs_lo = ALPHA_F*high(last,7) + ALPHA_C*high(last,8)
+        rho_max = U_SCRAM + ALPHA_F*(low(last,7) - T_F0) +
+                            ALPHA_C*(low(last,8) - T_C0)
+        rho_min = U_SCRAM + ALPHA_F*(high(last,7) - T_F0) +
+                            ALPHA_C*(high(last,8) - T_C0)
+        sdm_ok = sdm_lhs_hi <= SDM_RHS
+
         println("  n at T_probe (reconstructed): [$(round(n_final_lo; sigdigits=4)), $(round(n_final_hi; sigdigits=4))]")
         println("  T_c at T_probe: [$(round(Tc_final[1]; digits=2)), $(round(Tc_final[2]; digits=2))] °C")
         println("  T_f at T_probe: [$(round(Tf_final[1]; digits=2)), $(round(Tf_final[2]; digits=2))] °C")
+        println("  rho at T_probe: [$(round(rho_min; sigdigits=4)), $(round(rho_max; sigdigits=4))]  (shutdown margin = $(round(-rho_max; sigdigits=4)) dk/k)")
+        println("  shutdown_margin LHS sup: $(round(sdm_lhs_hi; sigdigits=4))  vs RHS $(round(SDM_RHS; sigdigits=4))  → $(sdm_ok ? "✓ DISCHARGED" : "× not yet")")
 
         results[T_probe] = (status="OK", n_sets=n_sets, elapsed=elapsed,
                             n_final=(n_final_lo, n_final_hi),
-                            Tc=Tc_final, Tf=Tf_final, Tcold=Tcold_final)
+                            Tc=Tc_final, Tf=Tf_final, Tcold=Tcold_final,
+                            sdm_lhs=(sdm_lhs_lo, sdm_lhs_hi),
+                            rho=(rho_min, rho_max),
+                            sdm_ok=sdm_ok,
+                            t_arr=t_arr,
+                            n_lo_ts=n_lo_ts,     n_hi_ts=n_hi_ts,
+                            rho_lo_ts=rho_lo_ts, rho_hi_ts=rho_hi_ts,
+                            Tc_lo_ts=Tc_lo_ts,   Tc_hi_ts=Tc_hi_ts,
+                            Tf_lo_ts=Tf_lo_ts,   Tf_hi_ts=Tf_hi_ts)
     catch err
         msg = sprint(showerror, err)
         println("  FAILED: ", first(msg, 300))
@@ -133,8 +191,8 @@ for T_probe in (10.0, 30.0, 60.0)
     haskey(results, T_probe) || continue
     r = results[T_probe]
     if r.status == "OK"
-        ok_subcrit = r.n_final[2] <= 1e-4 ? "✓ subcritical (n ≤ 1e-4)" : "× still above 1e-4"
+        ok_str = r.sdm_ok ? "✓ shutdown_margin DISCHARGED" : "× shutdown_margin not yet"
-        println("  T = $(T_probe) s: $(r.n_sets) sets, $(round(r.elapsed; digits=1))s wall — n ∈ [$(round(r.n_final[1]; sigdigits=3)), $(round(r.n_final[2]; sigdigits=3))]  $ok_subcrit")
+        println("  T = $(T_probe) s: $(r.n_sets) sets, $(round(r.elapsed; digits=1))s wall — rho ∈ [$(round(r.rho[1]; sigdigits=3)), $(round(r.rho[2]; sigdigits=3))]  $ok_str")
     else
         println("  T = $(T_probe) s: FAILED")
     end
@@ -154,7 +212,22 @@ for T_probe in (10.0, 30.0, 60.0)
         saved["T_$(Int(T_probe))_Tf_hi"]   = r.Tf[2]
         saved["T_$(Int(T_probe))_Tcold_lo"] = r.Tcold[1]
         saved["T_$(Int(T_probe))_Tcold_hi"] = r.Tcold[2]
+        saved["T_$(Int(T_probe))_sdm_lhs_hi"] = r.sdm_lhs[2]
+        saved["T_$(Int(T_probe))_rho_max"]    = r.rho[2]
+        saved["T_$(Int(T_probe))_sdm_ok"]     = r.sdm_ok ? 1.0 : 0.0
+        # Per-step time series for tube plotting.
+        saved["T_$(Int(T_probe))_t_arr"]     = r.t_arr
+        saved["T_$(Int(T_probe))_n_lo_ts"]   = r.n_lo_ts
+        saved["T_$(Int(T_probe))_n_hi_ts"]   = r.n_hi_ts
+        saved["T_$(Int(T_probe))_rho_lo_ts"] = r.rho_lo_ts
+        saved["T_$(Int(T_probe))_rho_hi_ts"] = r.rho_hi_ts
+        saved["T_$(Int(T_probe))_Tc_lo_ts"]  = r.Tc_lo_ts
+        saved["T_$(Int(T_probe))_Tc_hi_ts"]  = r.Tc_hi_ts
+        saved["T_$(Int(T_probe))_Tf_lo_ts"]  = r.Tf_lo_ts
+        saved["T_$(Int(T_probe))_Tf_hi_ts"]  = r.Tf_hi_ts
     end
 end
+saved["sdm_rhs"] = SDM_RHS
+saved["rho_sdm"] = RHO_SDM
 matwrite(mat_out, saved)
 println("\nSaved scram envelope summaries to $mat_out")

code/scripts/reach/reach_scram_pj_fat.jl

@@ -43,6 +43,10 @@ const T_F0    = T_C0 + P0 / HA
 const U_SCRAM    = -8 * BETA
 const Q_SG_DECAY = 0.03 * P0
 
+# X_exit threshold: shutdown_margin halfspace, mirrors predicates.json.
+const RHO_SDM = 0.01                                            # 1% dk/k
+const SDM_RHS = -RHO_SDM - U_SCRAM + ALPHA_F*T_F0 + ALPHA_C*T_C0   # ≈ 0.00297
+
 # Taylorized scram RHS (same as reach_scram_pj.jl).
 @taylorize function rhs_scram_fat_taylor!(dx, x, p, t)
     rho = U_SCRAM + ALPHA_F * (x[7] - T_F0) + ALPHA_C * (x[8] - T_C0)
@@ -198,6 +202,32 @@ for T_probe in (10.0, 30.0, 60.0)
         flow_hr = overapproximate(flow, Hyperrectangle)
         n_sets = length(flow_hr)
         println("  TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s")
+
+        # --- Per-step envelopes for plotting tubes ---
+        t_arr     = zeros(n_sets)
+        n_lo_ts   = zeros(n_sets); n_hi_ts   = zeros(n_sets)
+        rho_lo_ts = zeros(n_sets); rho_hi_ts = zeros(n_sets)
+        Tc_lo_ts  = zeros(n_sets); Tc_hi_ts  = zeros(n_sets)
+        Tf_lo_ts  = zeros(n_sets); Tf_hi_ts  = zeros(n_sets)
+        for (k, R) in enumerate(flow_hr)
+            s = set(R)
+            t_arr[k] = high(s, 10)
+            sumLC_lo_k = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
+                         LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
+            sumLC_hi_k = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
+                         LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
+            rho_lo_k = U_SCRAM + ALPHA_F*(high(s,7) - T_F0) + ALPHA_C*(high(s,8) - T_C0)
+            rho_hi_k = U_SCRAM + ALPHA_F*(low(s,7)  - T_F0) + ALPHA_C*(low(s,8)  - T_C0)
+            denom_lo_k = BETA - rho_hi_k
+            denom_hi_k = BETA - rho_lo_k
+            n_lo_ts[k]   = denom_lo_k > 0 ? LAMBDA * sumLC_lo_k / denom_hi_k : NaN
+            n_hi_ts[k]   = denom_lo_k > 0 ? LAMBDA * sumLC_hi_k / denom_lo_k : NaN
+            rho_lo_ts[k] = rho_lo_k
+            rho_hi_ts[k] = rho_hi_k
+            Tc_lo_ts[k]  = low(s, 8); Tc_hi_ts[k] = high(s, 8)
+            Tf_lo_ts[k]  = low(s, 7); Tf_hi_ts[k] = high(s, 7)
+        end
+
         last = set(flow_hr[end])
         sumLC_lo = LAM_1*low(last,1) + LAM_2*low(last,2) + LAM_3*low(last,3) +
                    LAM_4*low(last,4) + LAM_5*low(last,5) + LAM_6*low(last,6)
@@ -209,11 +239,25 @@ for T_probe in (10.0, 30.0, 60.0)
         denom_hi = BETA - rho_lo
         n_lo = denom_lo > 0 ? LAMBDA * sumLC_lo / denom_hi : NaN
         n_hi = denom_lo > 0 ? LAMBDA * sumLC_hi / denom_lo : NaN
+
+        # shutdown_margin discharge: alpha_f*T_f + alpha_c*T_c ≤ SDM_RHS.
+        # Coefficients negative → sup over the box at low(T_f), low(T_c).
+        sdm_lhs_hi = ALPHA_F*low(last,7) + ALPHA_C*low(last,8)
+        sdm_ok = sdm_lhs_hi <= SDM_RHS
+
         println("  n envelope:    [$(round(n_lo; sigdigits=4)), $(round(n_hi; sigdigits=4))]")
         println("  T_c envelope:  [$(round(low(last,8); digits=2)), $(round(high(last,8); digits=2))] °C")
         println("  T_f envelope:  [$(round(low(last,7); digits=2)), $(round(high(last,7); digits=2))] °C")
+        println("  rho envelope:  [$(round(rho_lo; sigdigits=4)), $(round(rho_hi; sigdigits=4))]  (shutdown margin = $(round(-rho_hi; sigdigits=4)) dk/k)")
+        println("  shutdown_margin LHS sup: $(round(sdm_lhs_hi; sigdigits=4))  vs RHS $(round(SDM_RHS; sigdigits=4))  → $(sdm_ok ? "✓ DISCHARGED" : "× not yet")")
 
-        results[T_probe] = (status="OK", n=(n_lo, n_hi), elapsed=elapsed)
+        results[T_probe] = (status="OK", n=(n_lo, n_hi), elapsed=elapsed,
+                            rho=(rho_lo, rho_hi), sdm_lhs_hi=sdm_lhs_hi, sdm_ok=sdm_ok,
+                            t_arr=t_arr,
+                            n_lo_ts=n_lo_ts,     n_hi_ts=n_hi_ts,
+                            rho_lo_ts=rho_lo_ts, rho_hi_ts=rho_hi_ts,
+                            Tc_lo_ts=Tc_lo_ts,   Tc_hi_ts=Tc_hi_ts,
+                            Tf_lo_ts=Tf_lo_ts,   Tf_hi_ts=Tf_hi_ts)
     catch err
         println("  FAILED: ", first(sprint(showerror, err), 300))
         results[T_probe] = (status="FAILED",)
@@ -226,7 +270,8 @@ for T_probe in (10.0, 30.0, 60.0)
     haskey(results, T_probe) || continue
     r = results[T_probe]
     if r.status == "OK"
-        @printf "  T = %4.0f s: n ∈ [%.3e, %.3e]\n" T_probe r.n[1] r.n[2]
+        ok_str = r.sdm_ok ? "✓ shutdown_margin DISCHARGED" : "× shutdown_margin not yet"
+        @printf "  T = %4.0f s: n ∈ [%.3e, %.3e]  ρ ∈ [%.4f, %.4f]  %s\n" T_probe r.n[1] r.n[2] r.rho[1] r.rho[2] ok_str
     else
         println("  T = $T_probe s: FAILED")
     end
@@ -234,11 +279,26 @@ end
 
 mat_out = joinpath(results_dir, "reach_scram_pj_fat.mat")
 saved = Dict{String,Any}("fat_lo" => fat_lo, "fat_hi" => fat_hi,
-                         "sources" => ["shutdown", "heatup_tight", "operation", "loca_operation"])
+                         "sources" => ["shutdown", "heatup_tight", "operation", "loca_operation"],
+                         "sdm_rhs" => SDM_RHS, "rho_sdm" => RHO_SDM)
 for (T_probe, r) in results
     if r.status == "OK"
         saved["T_$(Int(T_probe))_n_lo"] = r.n[1]
         saved["T_$(Int(T_probe))_n_hi"] = r.n[2]
+        saved["T_$(Int(T_probe))_rho_lo"] = r.rho[1]
+        saved["T_$(Int(T_probe))_rho_hi"] = r.rho[2]
+        saved["T_$(Int(T_probe))_sdm_lhs_hi"] = r.sdm_lhs_hi
+        saved["T_$(Int(T_probe))_sdm_ok"] = r.sdm_ok ? 1.0 : 0.0
+        # Per-step time series for tube plotting.
+        saved["T_$(Int(T_probe))_t_arr"]     = r.t_arr
+        saved["T_$(Int(T_probe))_n_lo_ts"]   = r.n_lo_ts
+        saved["T_$(Int(T_probe))_n_hi_ts"]   = r.n_hi_ts
+        saved["T_$(Int(T_probe))_rho_lo_ts"] = r.rho_lo_ts
+        saved["T_$(Int(T_probe))_rho_hi_ts"] = r.rho_hi_ts
+        saved["T_$(Int(T_probe))_Tc_lo_ts"]  = r.Tc_lo_ts
+        saved["T_$(Int(T_probe))_Tc_hi_ts"]  = r.Tc_hi_ts
+        saved["T_$(Int(T_probe))_Tf_lo_ts"]  = r.Tf_lo_ts
+        saved["T_$(Int(T_probe))_Tf_hi_ts"]  = r.Tf_hi_ts
     end
 end
 matwrite(mat_out, saved)

code/scripts/sim/main_mode_sweep.jl

@@ -17,32 +17,23 @@ using MatrixEquations
 using JSON
 
 include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_states.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_solver.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "plot_pke_results.jl"))
 include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
 
-# --- Plant + predicates ---
+# --- Plant + predicates + canonical ICs ---
 plant = pke_params()
 
-# Load T_standby from predicates.json for the hot-standby IC + heatup ref.
 pred_path = joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")
-pred_raw  = JSON.parsefile(pred_path)
+predicates = JSON.parsefile(pred_path)
-T_standby = plant.T_c0 + pred_raw["derived"]["T_standby_offset_C"]
+states = pke_states(plant; predicates=predicates)
-
+T_standby = states.T_standby
-# --- ICs ---
+x_op   = states.op
-x_op = pke_initial_conditions(plant)
+x_shut = states.shutdown
-
+x_heat = states.heatup
-# Hot-standby IC: everything at T_standby, trace n.
-n_shut = 1e-6
-C_shut = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_shut
-x_shut = [n_shut; C_shut; T_standby; T_standby; T_standby]
-
-# Heatup IC: already critical at 0.1% power, thermally matched.
-n_heat = 1e-3
-C_heat = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_heat
-x_heat = [n_heat; C_heat; T_standby; T_standby; T_standby]
 
 # --- LQR gain (needed by ctrl_operation_lqr_factory) ---
 A, B, B_w, _, _, _ = pke_linearize(plant)

code/scripts/sim/validate_pj.jl

@@ -19,22 +19,22 @@ using OrdinaryDiffEq
 using Plots
 
 include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
+include(joinpath(@__DIR__, "..", "..", "src", "pke_states.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
 include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs_pj.jl"))
 include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
 
 plant = pke_params()
-T_standby = plant.T_c0 - 33.333333
+states = pke_states(plant)   # falls back to default offset; no predicates needed for validation
+T_standby = states.T_standby
 
 # Heatup scenario — same as sim_heatup.jl.
 ref_heat = (; T_start=T_standby, T_target=plant.T_c0, ramp_rate=28/3600)
 Q_sg = t -> 0.0
 
-# Initial conditions — both at n=1e-3, same temperatures.
+# Initial conditions — both at n=1e-3 (the heatup canonical IC), same temperatures.
-n0 = 1e-3
+x0_full = states.heatup            # 10-state
-C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n0
+x0_pj   = states.heatup[2:end]     # drop n; PJ reconstructs it
-x0_full = [n0; C0; T_standby; T_standby; T_standby]
-x0_pj   = [C0; T_standby; T_standby; T_standby]
 
 tspan = (0.0, 3000.0)
 

code/src/pke_states.jl

@@ -0,0 +1,61 @@
+"""
+    pke_states(plant; predicates=nothing,
+               T_standby_offset_C=-33.333333,
+               n_shutdown=1e-6, n_heatup=1e-3)
+
+Single-point initial conditions for every canonical DRC mode.
+
+Returns a NamedTuple with:
+
+- `T_standby` — scalar, hot-standby coolant temperature [°C]
+- `op`       — full-power steady state (`n = 1`, precursors at equilibrium)
+- `shutdown` — hot standby (trace `n`, all temperatures at `T_standby`)
+- `heatup`   — critical at low power, thermals matched to `T_standby`
+
+Each state is a 10-element `Vector{Float64}` matching
+`[n; C₁..C₆; T_f; T_c; T_cold]`.
+
+# Source of truth for `T_standby`
+
+Pass `predicates=JSON.parsefile("reachability/predicates.json")` to read
+the canonical offset from `predicates.derived.T_standby_offset_C`. If
+omitted, falls back to the keyword default (`-33.333333` °C, equivalent
+to `-60` °F — matches the `predicates.json` value at the time of writing).
+
+# Why a single function
+
+ICs were previously duplicated across `main_mode_sweep.jl`,
+`barrier_sos_2d_shutdown.jl`, `validate_pj.jl`, etc. Each copy independently
+defined `T_standby`, `n_shut = 1e-6`, `n_heat = 1e-3`, and the precursor-
+equilibrium formula. Lifting once means you change power-level conventions
+in one place.
+
+# X_entry polytopes
+
+This function returns single-point ICs for forward simulation and
+equilibrium finding. The `X_entry` *polytopes* used by reach analysis
+live in `reachability/predicates.json::mode_definitions.<mode>.X_entry_polytope`
+and should be loaded directly by reach scripts — don't conflate the two.
+"""
+function pke_states(plant; predicates=nothing,
+                    T_standby_offset_C=-33.333333,
+                    n_shutdown=1e-6, n_heatup=1e-3)
+    if predicates !== nothing
+        T_standby = plant.T_c0 + predicates["derived"]["T_standby_offset_C"]
+    else
+        T_standby = plant.T_c0 + T_standby_offset_C
+    end
+
+    # Full-power steady state.
+    op = pke_initial_conditions(plant)
+
+    # Hot-standby: trace n, all thermals at T_standby.
+    C_shut = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_shutdown
+    shutdown = [n_shutdown; C_shut; T_standby; T_standby; T_standby]
+
+    # Heatup: critical at 0.1% power, thermals matched to T_standby.
+    C_heat = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_heatup
+    heatup = [n_heatup; C_heat; T_standby; T_standby; T_standby]
+
+    return (; T_standby, op, shutdown, heatup)
+end

code/src/sos_barrier.jl

@@ -0,0 +1,118 @@
+"""
+    solve_sos_barrier_2d(A_red, x_vars, entry_halfspaces, unsafe_halfspaces;
+                         barrier_degree=4, multiplier_degree=2,
+                         eps_cap=1.0, solver_attrs=("printlevel" => 0))
+
+Find a degree-`barrier_degree` polynomial barrier certificate `B(x)` for
+the linear closed-loop `dx/dt = A_red x` on a 2-D state projection,
+using the Prajna–Jadbabaie SOS formulation.
+
+A *barrier certificate* `B(x)` proves forward invariance of `{B(x) ≤ 0}`:
+- `B ≤ -ε` on `X_entry` (the certified-safe entry set)
+- `B ≥ +ε` on each unsafe halfspace
+- `dB/dt = ∇B · A_red·x ≤ 0` everywhere
+
+If such a `B` exists with `ε > 0`, the entry set never reaches the
+unsafe set under the closed-loop flow.
+
+# Why ε is here
+
+The naive feasibility version (no slack, no objective) silently returns
+`B ≡ 0` — vacuously satisfies `B ≤ 0` and `B ≥ 0` with `σ = 0`. We
+maximize `ε ∈ [0, eps_cap]` to force strict separation; `eps_cap` is
+finite because `B` has free scale (without a cap the primal is
+unbounded → `DUAL_INFEASIBLE`). We only care whether `ε* > 0`; its
+magnitude is just the unit of `B`.
+
+# Arguments
+
+- `A_red::AbstractMatrix` — 2×2 closed-loop dynamics matrix.
+- `x_vars` — tuple `(x1, x2)` of `DynamicPolynomials.PolyVar`s. Caller
+  must `@polyvar` these so the polynomial expressions in
+  `entry_halfspaces` and `unsafe_halfspaces` reference the same variables.
+- `entry_halfspaces::Vector{<:AbstractPolynomialLike}` — list of
+  polynomials `g_e(x)` such that `X_entry = ⋂ {g_e(x) ≥ 0}`.
+- `unsafe_halfspaces::Vector{<:AbstractPolynomialLike}` — list of
+  polynomials `g_u(x)` such that each `X_unsafe_i = {g_u_i(x) ≥ 0}`.
+  The certificate proves separation from the *union* of these regions.
+
+# Keyword arguments
+
+- `barrier_degree::Int=4` — degree of the barrier polynomial `B`.
+- `multiplier_degree::Int=2` — degree of the SOS multiplier polynomials.
+- `eps_cap::Real=1.0` — upper bound on `ε` (just controls the scale of
+  the returned `B`).
+- `solver_attrs` — passed to `optimizer_with_attributes(CSDP.Optimizer, ...)`.
+
+# Returns
+
+NamedTuple `(; status, B, ε, model)`:
+- `status` — `MOI.TerminationStatusCode`. `OPTIMAL` with `ε > 1e-8` means
+  a real certificate.
+- `B` — the optimized polynomial value, or `nothing` on failure.
+- `ε` — the achieved separation slack, or `nothing` on failure.
+- `model` — the JuMP `SOSModel`, returned for callers that want to
+  inspect multipliers.
+
+# Caveats
+
+- The Lie-derivative condition is checked *globally* (∇B·f ∈ SOS), not
+  just on `{B = 0}`. This is convex but stronger than needed; it
+  effectively requires `B` to be a Lyapunov function. Use the bilinear
+  Putinar form for a tighter test (out of scope here — needs alternation).
+- No disturbance term. Extend by adding `B_w·w` worst-case to the Lie
+  inequality and using a Putinar multiplier on the disturbance polytope.
+- The 2-D projection is a modeling simplification — cross-coupling from
+  dropped states is not certified. Norm-check the dropped block at the
+  caller site.
+"""
+function solve_sos_barrier_2d(A_red, x_vars, entry_halfspaces, unsafe_halfspaces;
+                              barrier_degree::Int=4,
+                              multiplier_degree::Int=2,
+                              eps_cap::Real=1.0,
+                              solver_attrs=("printlevel" => 0,))
+    x1, x2 = x_vars
+    f_nom = A_red * [x1; x2]
+
+    solver = optimizer_with_attributes(CSDP.Optimizer, solver_attrs...)
+    model = SOSModel(solver)
+
+    monos_B = monomials([x1, x2], 0:barrier_degree)
+    @variable(model, B_poly, Poly(monos_B))
+
+    monos_σ = monomials([x1, x2], 0:multiplier_degree)
+
+    # Slack variable: forces strict separation. Cap is finite because B
+    # has free scale (otherwise primal unbounded → DUAL_INFEASIBLE).
+    @variable(model, 0 <= ε <= eps_cap)
+    @objective(model, Max, ε)
+
+    # (a) B + ε ≤ 0 on X_entry.
+    σ_e_terms = sum(begin
+        σ = @variable(model, [1:1], SOSPoly(monos_σ))[1]
+        σ * g
+    end for g in entry_halfspaces)
+    @constraint(model, -B_poly - ε - σ_e_terms in SOSCone())
+
+    # (b) B - ε ≥ 0 on each X_unsafe_i (separately — union semantics).
+    for g_u in unsafe_halfspaces
+        σ_u = @variable(model, [1:1], SOSPoly(monos_σ))[1]
+        @constraint(model, B_poly - ε - σ_u * g_u in SOSCone())
+    end
+
+    # (c) Lie derivative ≤ 0 globally.
+    dB_dx = [differentiate(B_poly, x1), differentiate(B_poly, x2)]
+    lie = dB_dx[1] * f_nom[1] + dB_dx[2] * f_nom[2]
+    @constraint(model, -lie in SOSCone())
+
+    optimize!(model)
+    status = termination_status(model)
+
+    if status == MOI.OPTIMAL
+        ε_val = value(ε)
+        B_val = value(B_poly)
+        return (; status, B=B_val, ε=ε_val, model)
+    else
+        return (; status, B=nothing, ε=nothing, model)
+    end
+end

docs/model_cheatsheet.md

@@ -0,0 +1,220 @@
+# PWR Plant Model Cheatsheet
+
+Quick reference for the 10-state PKE + 3-node thermal-hydraulics model
+used throughout this demo.
+
+Source of truth: `code/src/pke_params.jl` (constants + steady state),
+`code/src/pke_th_rhs.jl` (dynamics), `code/src/pke_th_rhs_pj.jl`
+(prompt-jump variant), `code/controllers/controllers.jl` (control laws).
+
+All values in SI internally; printed/plotted in °F.
+
+---
+
+## State vector
+
+`x = [n, C₁, C₂, C₃, C₄, C₅, C₆, T_f, T_c, T_cold]`  (10 states)
+
+| Symbol | Index | Meaning | Units |
+|---|---|---|---|
+| n | 1 | Normalized neutron population (n=1 ⇔ full power) | — |
+| Cᵢ | 2–7 | Delayed-neutron precursor concentrations (6 groups) | — |
+| T_f | 8 | Fuel-pellet bulk temperature | °C |
+| T_c | 9 | Average coolant temperature in core (= T_avg) | °C |
+| T_cold | 10 | Cold-leg coolant temperature | °C |
+
+**Derived (not states):**
+- `T_hot = 2·T_c − T_cold` (linear coolant profile assumption)
+- `P = P₀·n` (thermal power)
+
+**Inputs:**
+- `u` = control input (rod-induced reactivity), units of Δk/k
+- `Q_sg(t)` = bounded disturbance (steam-generator heat removal), W
+
+---
+
+## Reactivity (the coupling between neutronics and thermal)
+
+$$\rho(x, u) \;=\; u \;+\; \alpha_f\,(T_f - T_{f,0}) \;+\; \alpha_c\,(T_c - T_{c,0})$$
+
+Three contributions:
+- `u` — externally controlled (rod motion)
+- `α_f·(T_f − T_f0)` — Doppler feedback (negative; fuel heats → more U-238 absorption → less reactivity)
+- `α_c·(T_c − T_c0)` — moderator-density feedback (negative; coolant heats → less moderation → less reactivity)
+
+Both feedback coefficients are negative — the plant is **self-stabilizing**.
+
+---
+
+## ODE system (full 10-state)
+
+### Neutronics (Point-Kinetic Equations)
+
+$$\frac{dn}{dt} \;=\; \frac{\rho - \beta}{\Lambda}\,n \;+\; \sum_{i=1}^{6} \lambda_i C_i$$
+
+$$\frac{dC_i}{dt} \;=\; \frac{\beta_i}{\Lambda}\,n \;-\; \lambda_i C_i \qquad (i = 1, \ldots, 6)$$
+
+The first equation has a **stiff coefficient** `1/Λ ≈ 10⁴ s⁻¹` — this is what makes the full plant unverifiable by direct nonlinear reach.
+
+### Thermal-hydraulics (3 lumped nodes)
+
+$$\frac{dT_f}{dt} \;=\; \frac{P_0\,n \;-\; hA\,(T_f - T_c)}{M_f\,c_f}$$
+
+$$\frac{dT_c}{dt} \;=\; \frac{hA\,(T_f - T_c) \;-\; 2\,W\,c_c\,(T_c - T_{\text{cold}})}{M_c\,c_c}$$
+
+$$\frac{dT_{\text{cold}}}{dt} \;=\; \frac{2\,W\,c_c\,(T_c - T_{\text{cold}}) \;-\; Q_{sg}}{M_{sg}\,c_c}$$
+
+Physical chain: `fission heat → fuel → coolant → cold leg → SG → back to cold leg`.
+
+The factor of 2 in `2·W·c_c·(T_c − T_cold)` comes from the linear-profile
+assumption: the enthalpy carried by flow at avg `T_c` from inlet `T_cold`
+to outlet `T_hot = 2T_c − T_cold` is `W·c_c·(T_hot − T_cold) = 2·W·c_c·(T_c − T_cold)`.
+
+---
+
+## Prompt-jump (PJ) reduction
+
+Set `dn/dt = 0` and solve algebraically for `n`:
+
+$$n(x) \;=\; \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho(x, u)}$$
+
+State drops 10 → 9 (drop `n`; reconstruct it from C and ρ). Stiffness gone.
+Validity condition: `β − ρ > 0` with margin (encoded as the
+`prompt_critical_margin_*` halfspace per mode).
+
+Tikhonov bound: `|x(t) − x_PJ(t)| ≤ C·Λ = O(10⁻⁴)` — characterized
+residual error vs the 10-state plant.
+
+---
+
+## Constants
+
+### Neutronics (6-group point-kinetic)
+
+| Symbol | Value | Meaning |
+|---|---|---|
+| Λ | `1×10⁻⁴ s` | Prompt-neutron generation time |
+| β = Σβᵢ | `0.006502` | Total delayed-neutron fraction (sum of 6 groups) |
+| β₁..β₆ | `[0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273]` | Per-group delayed fractions |
+| λ₁..λ₆ | `[0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01]` s⁻¹ | Precursor decay constants |
+
+(Standard 6-group U-235 thermal-fission values from Keepin.)
+
+### Thermal-hydraulic
+
+| Symbol | Value | Meaning |
+|---|---|---|
+| P₀ | `1×10⁹ W` | Rated thermal power (1 GWth) |
+| M_f | `50 000 kg` | Fuel mass (lumped) |
+| c_f | `300 J/(kg·K)` | Fuel specific heat |
+| M_c | `20 000 kg` | Core coolant mass |
+| c_c | `5 450 J/(kg·K)` | Coolant specific heat |
+| hA | `5×10⁷ W/K` | Fuel-to-coolant heat-transfer coefficient |
+| W | `5 000 kg/s` | Primary coolant mass flow rate |
+| M_sg | `30 000 kg` | Steam-generator coolant inventory |
+
+### Reactivity feedback coefficients (linearized about full-power op point)
+
+| Symbol | Value | Meaning |
+|---|---|---|
+| α_f | `−2.5×10⁻⁵ /°C` | Fuel/Doppler temperature coefficient |
+| α_c | `−1×10⁻⁴ /°C` | Moderator/coolant temperature coefficient |
+
+These are linear slopes about the full-power operating point. **Trust
+range: ~±50 °C around T_f0/T_c0.** Cold-shutdown extrapolation breaks.
+
+### Derived steady-state (full-power equilibrium)
+
+Single free parameter: `T_cold0 = 290 °C` (full-power cold-leg temp).
+Everything else falls out of energy balance:
+
+| Quantity | Derivation | Value |
+|---|---|---|
+| ΔT_core | `P₀ / (W·c_c)` | `36.7 °C` |
+| T_hot0 | `T_cold0 + ΔT_core` | `326.7 °C` |
+| T_c0 (= T_avg0) | `(T_hot0 + T_cold0) / 2` | `308.35 °C` |
+| T_f0 | `T_c0 + P₀/hA` | `328.35 °C` |
+| n₀ | (definition) | `1.0` |
+| C_i0 | `βᵢ·n₀ / (λᵢ·Λ)` | (from `dC/dt = 0`) |
+
+### Hot-standby reference
+
+| Quantity | Derivation | Value |
+|---|---|---|
+| T_standby | `T_c0 − 33.33 °C` (`= T_c0 − 60 °F`) | `275.02 °C` |
+
+This is the canonical "hot standby" operating point — coolant warm but
+power below criticality. Defined as a 60 °F offset below full-power T_avg.
+
+---
+
+## Per-mode controllers
+
+| Mode | Control law | Constants | Notes |
+|---|---|---|---|
+| **`q_shutdown`** | `u = −5β` (constant rod) | — | Open-loop rod-held |
+| **`q_heatup`** | `u = K_p·(T_ref(t) − T_c)` (FL/proportional) | `K_p = 1×10⁻⁴` | T_ref ramps from T_standby at 28 °C/hr (tech-spec heatup limit) |
+| **`q_operation` (P)** | `u = K_p·(T_avg_ref − T_avg)` | `K_p ≈ 1×10⁻⁴` | Simple proportional T_avg tracker |
+| **`q_operation` (LQR)** | `u = −K_LQR·δx` | gain solved from Riccati | Cached in factory closure |
+| **`q_scram`** | `u = −8β` (constant rod, max insertion) | — | Open-loop rod-held; full negative reactivity |
+
+Heatup ramp rate: `dT_ref/dt = 28/3600 ≈ 0.00778 °C/s` (28 °C/hr — US PWR tech-spec maximum heatup rate, set by vessel thermal-stress limits).
+
+---
+
+## Key predicates / halfspaces (used by reach analysis)
+
+From `reachability/predicates.json`:
+
+| Predicate | Concretization | What it discharges |
+|---|---|---|
+| `t_avg_above_min` | `T_c ≥ T_standby + 5.556 °C` (= 280.58 °C) | Shutdown→heatup transition trigger |
+| `t_avg_in_range` | `\|T_c − T_c0\| ≤ 2.78 °C` (= [305.57, 311.13]) | Heatup→operation transition trigger |
+| `p_above_crit` | `n ≥ 1×10⁻⁴` | "Reactor at-or-above critical" (also part of heatup→operation) |
+| `fuel_centerline` | `T_f ≤ 1200 °C` | UO₂ melt prevention |
+| `t_avg_high_trip` | `T_c ≤ 320 °C` | Reactor trip limit |
+| `t_avg_low_trip` | `T_c ≥ 280 °C` | Reactor trip limit |
+| `n_high_trip` | `n ≤ 1.15` | High-flux trip (118% nominal) |
+| `cold_leg_subcooled` | `T_cold ≤ T_cold0 + 15 = 305 °C` | Subcooling margin |
+| `heatup_rate_upper` | `0.4587·T_f − 0.9587·T_c + 0.5·T_cold ≤ 0.01389 °C/s` | Coolant heatup ≤ 50 °C/hr (28 + overshoot) |
+| `heatup_rate_lower` | (mirror, lower bound) | Cooldown ≤ 50 °C/hr |
+| `prompt_critical_margin_*` | `β − ρ ≥ δ` (per-mode form) | PJ reduction validity |
+| `shutdown_margin` | `α_f·T_f + α_c·T_c ≤ 0.00297` | Scram success: rho ≤ −0.01 (1% Δk/k) |
+
+The `dT_c/dt` halfspace coefficients above come from differentiating
+the `T_c` ODE: `a_f = hA/(M_c·c_c) = 0.4587 /s`,
+`a_c = −(hA + 2W·c_c)/(M_c·c_c) = −0.9587 /s`,
+`a_cold = 2W·c_c/(M_c·c_c) = 0.5 /s`. Sums to zero by equilibrium.
+
+---
+
+## Per-mode obligations
+
+| Mode | Kind | Obligation |
+|---|---|---|
+| `q_shutdown` | equilibrium | Stay in X_safe forever; transition out when `t_avg_above_min` becomes true |
+| `q_heatup` | transition | From X_entry, reach X_exit within `[T_min=7714, T_max=18000]` s, maintain `inv1_holds` |
+| `q_operation` | equilibrium | Stay in X_safe forever under bounded `Q_sg` |
+| `q_scram` | transition | From any state, drive to `shutdown_margin` within `T_max=60` s, maintain bounded T |
+
+`inv1_holds` (heatup safety invariant) = conjunction of `fuel_centerline`,
+`cold_leg_subcooled`, `heatup_rate_upper`, `heatup_rate_lower`,
+`prompt_critical_margin_heatup`.
+
+---
+
+## Sanity numbers (to memorize before the talk)
+
+| Quantity | Value |
+|---|---|
+| Full-power n | `1.0` |
+| Full-power T_avg | `308 °C ≈ 587 °F` |
+| Hot-standby T_avg | `275 °C ≈ 527 °F` |
+| Heatup span (standby → operation) | `33 °C ≈ 60 °F` |
+| Tech-spec heatup rate | `28 °C/hr` |
+| Nominal heatup time | `33 / 28 = 1.19 hr ≈ 71 min` |
+| T_min (heatup obligation) | `7 714 s = 2 hr 8 min` |
+| T_max (heatup obligation) | `18 000 s = 5 hr` |
+| Stiffness ratio (full plant) | `~10⁵` (Λ⁻¹ ÷ thermal time-constant) |
+| Scram rod worth | `−8β = −0.052` Δk/k (~5% subcritical) |
+| Shutdown-margin requirement | `\|ρ\| ≥ 0.01 = 1% Δk/k` |

journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex

@@ -0,0 +1,217 @@
+% ---------------------------------------------------------------------------
+% 2026-04-27 — Hot-standby SOS barrier + scram X_exit redefinition
+% Live / B-style entry: two task closeouts in one sitting.
+% ---------------------------------------------------------------------------
+
+\session{2026-04-27}{Hacker-Split, Mac mini}{%
+Two thesis-relevant closeouts in one go: redefine the scram exit
+predicate from a power threshold to a shutdown-margin halfspace
+(actual NRC criterion, and a clean linear halfspace), then port the
+degree-4 SOS barrier machinery from operation to hot-standby and
+discharge a forever-invariance certificate on a 2-D thermal projection.}
+
+\section{2026-04-27 --- Scram exit \& hot-standby SOS}
+\label{sec:20260427-shutdown-sos-scram-exit}
+
+\subsection*{Scram \texttt{X\_exit}: from $n \leq 10^{-4}$ to shutdown margin}
+
+The scram-mode exit predicate as written in
+\texttt{reachability/predicates.json} read
+
+\[
+  X_{\mathrm{exit}}^{(\mathrm{scram})} \;=\; \{\, x \;:\; n \leq 10^{-4} \;\wedge\; T_f \leq T_{f,0} + 50\,^\circ\mathrm{C}\,\}.
+\]
+
+Two structural problems with this:
+
+\begin{enumerate}
+  \item \textbf{$n$ is nonlinear in the prompt-jump (PJ) state.} Under PJ,
+    $n = \Lambda \sum_i \lambda_i C_i / (\beta - \rho)$, and $\rho$
+    depends on $T_f, T_c$. So $n \leq 10^{-4}$ is not a halfspace in the
+    9-state PJ vector; the existing \texttt{reach\_scram\_pj.jl} was
+    reconstructing $n$ post-hoc and checking the bound on the endpoint
+    only.
+  \item \textbf{The $T_f$ bound is infeasible by 60\,s under decay heat.}
+    With $Q_{\mathrm{sg}} = 0.03 P_0$ and a fuel time constant
+    $M_f c_f / hA \approx 0.3\,\mathrm{s}$, the fuel rapidly equilibrates
+    with $T_c$ and the system loses only $\sim 5\,^\circ\mathrm{C}$ in
+    $60\,\mathrm{s}$. The conjoined bound was decorative.
+\end{enumerate}
+
+The actual NRC tech-spec criterion for scram success is phrased in
+shutdown margin ($\rho \leq -\rho_{\mathrm{SDM}}$, typically $1\%
+\Delta k/k$). Under constant rod $u = U_{\mathrm{SCRAM}} = -8\beta$,
+total reactivity is
+
+\[
+  \rho(x) \;=\; U_{\mathrm{SCRAM}} \;+\; \alpha_f (T_f - T_{f,0}) \;+\; \alpha_c (T_c - T_{c,0}),
+\]
+
+\emph{linear} in $(T_f, T_c)$ — a single-row halfspace in the reach state.
+Replacing $X_{\mathrm{exit}}$ with the predicate
+\texttt{shutdown\_margin}:
+
+\[
+  \alpha_f T_f + \alpha_c T_c \;\leq\;
+    -\rho_{\mathrm{SDM}} - U_{\mathrm{SCRAM}} + \alpha_f T_{f,0} + \alpha_c T_{c,0}
+    \;\approx\; 2.97 \times 10^{-3}.
+\]
+
+Re-ran \texttt{reach\_scram\_pj.jl} (TMJets, orderT $= 4$, orderQ $= 2$).
+Discharged at all three probe horizons:
+
+\begin{center}
+\begin{tabular}{r r r r r}
+$T$ (s) & reach-sets & wall (s) & $\rho$ at $T$ & discharged \\
+\hline
+$10$ & $6919$ & $98.6$ & $[-0.0507,\,-0.0504]$ & \checkmark \\
+$30$ & $9900$ & $130.5$ & $[-0.0506,\,-0.0503]$ & \checkmark \\
+$60$ & $12340$ & $164.2$ & $[-0.0503,\,-0.0500]$ & \checkmark \\
+\end{tabular}
+\end{center}
+
+Required $|\rho| \geq 0.01$. Actual $|\rho| \approx 0.05$ — five times
+the requirement, dominated by rod worth (as expected; Doppler/moderator
+contributions vary by only $\sim 3\%$ of $U_{\mathrm{SCRAM}}$).
+
+\begin{decision}
+\textbf{Shutdown margin (in $\Delta k/k$) is the canonical scram success
+criterion across all modes.} The $n$-threshold framing was a
+back-translation that didn't survive scrutiny. Where reach scripts
+need to discharge a successor obligation post-scram, use
+\texttt{shutdown\_margin}.
+\end{decision}
+
+Stale-constant gotcha caught while drafting the predicate: I derived
+the RHS by hand using rounded $T_{f,0} = 320$, $T_{c,0} = 300$, but
+the actual values from \texttt{pke\_params} are $T_{c,0} = 308.35$
+($= (T_{\mathrm{hot},0} + T_{\mathrm{cold},0})/2$ with
+$T_{\mathrm{hot},0} = T_{\mathrm{cold},0} + P_0 / (W c_c) = 326.7$)
+and $T_{f,0} = 328.35$. Off by $\Delta\mathrm{RHS} \sim 10^{-3}$. Fixed
+by switching the JSON \texttt{rhs\_expr} to a symbolic form that gets
+substituted at load time. Lesson: every constant in
+\texttt{predicates.json} that's derivable from \texttt{pke\_params}
+should be symbolic, not baked.
+
+\subsection*{SOS barrier on hot-standby (q\_shutdown)}
+
+Companion to \texttt{barrier\_sos\_2d.jl} (operation/LQR). New script
+\texttt{barrier\_sos\_2d\_shutdown.jl}. Hot-standby is an
+\emph{equilibrium-mode} obligation (forever-invariance), not
+reach-avoid, so SOS suits it well.
+
+Controller: $u_{\mathrm{shutdown}} = -5\beta = -0.0325$, constant. With
+$Q_{\mathrm{sg}} = 0$ (no SG load), the thermal subsystem is adiabatic:
+\emph{any} uniform temperature $T_f = T_c = T_{\mathrm{cold}} = T^*$
+sufficient to make $\rho < 0$ is invariant. The hot-standby IC sets
+$T^* = T_{\mathrm{standby}} = 275.02\,^\circ\mathrm{C}$; a 50\,000\,s sim
+confirms the trajectory parks there, $\|dx/dt\| \sim 10^{-22}$. So the
+linearization is at $x_{\mathrm{eq}} = (T_{\mathrm{standby}}, T_{\mathrm{standby}}, T_{\mathrm{standby}})$
+with $u_{\mathrm{eq}} = -5\beta$, $Q_* = 0$.
+
+\textbf{2-D reduction.} Picked $(T_c, T_{\mathrm{cold}})$ — the slow safety-relevant
+thermal modes. $n$ and the precursors are decoupled at quasi-zero power
+and not the safety driver in this mode. $T_f$ tracks $T_c$ with time
+constant $\sim 0.3\,\mathrm{s}$ and is folded into the dynamics implicitly
+(the 2-D reduction loses $\|cross\| = 0.459$ of coupling from dropped
+states — non-trivial; see \emph{Open issues} below).
+
+Reduced closed-loop:
+\[
+  A_{\mathrm{red}} \;=\; \begin{bmatrix} -0.959 & 0.500 \\ 0.333 & -0.333 \end{bmatrix},
+  \quad \mathrm{eig}(A_{\mathrm{red}}) = \{-1.16,\, -0.132\}.
+\]
+
+Both Hurwitz; slow mode is $T_{\mathrm{cold}}$ ($\tau \sim 7.6\,\mathrm{s}$).
+
+\textbf{Sets.} Safety: $|\delta T_c| \leq 10$, $|\delta T_{\mathrm{cold}}| \leq 15$.
+Entry: $|\delta T_c| \leq 5$, $|\delta T_{\mathrm{cold}}| \leq 5$ (matches the
+\texttt{q\_shutdown.X\_entry\_polytope} extent recentered on $x_{\mathrm{eq}}$).
+Unsafe focus: $\delta T_c \geq +10$ — over-warming is the harder direction
+because the controller has rods already maxed in negative reactivity and
+no recourse but to wait for the moderator coefficient.
+
+\textbf{Methodological gotcha: trivial barrier.} First run with the
+Prajna--Jadbabaie SOS feasibility formulation returned
+$B(x) \equiv 0$ — vacuously satisfies $B \leq 0$ on entry, $B \geq 0$
+on unsafe (with $\sigma_u \equiv 0$), $\nabla B \cdot f = 0$ globally.
+Standard fix: add a \emph{strict-separation slack} $\varepsilon > 0$
+and tighten the constraints to $B \leq -\varepsilon$ on entry,
+$B \geq +\varepsilon$ on unsafe. Maximize $\varepsilon$.
+
+Second run: \texttt{DUAL\_INFEASIBLE}. Because $B$ has free scale, the
+program is unbounded — scaling $B \to cB$ scales $\varepsilon \to c\varepsilon$.
+Cap $\varepsilon \leq 1$ (the unit is arbitrary; we only need
+$\varepsilon^* > 0$ for a real certificate).
+
+Third run: \texttt{OPTIMAL}, $\varepsilon^* = 1.0$ (hit the cap, meaning
+the primal is feasible with arbitrarily large separation modulo scale).
+Dropping numerical-noise terms,
+
+\[
+  B(x_1, x_2) \approx -16.91 \;+\; 0.022\,x_1^2 \;+\; 0.027\,x_2^2 \;+\; 0.005\,x_1^4 \;+\; \text{(cross / cubic terms)},
+\]
+
+with $x_1 = \delta T_c$, $x_2 = \delta T_{\mathrm{cold}}$. Geometry: a bowl
+with floor $B(0,0) = -16.91 \leq -\varepsilon$, rising past zero around
+$|x_1| \approx 7.5$, comfortably $\geq +\varepsilon$ at the unsafe
+threshold $x_1 = 10$ (where $B \approx +35$). $\dot B \leq 0$ globally
+because $A_{\mathrm{red}}$ is Hurwitz; the polynomial barrier is
+essentially a degree-4 Lyapunov function for this Hurwitz system.
+
+\begin{decision}
+SOS feasibility programs without an objective + scale normalization
+silently return $B \equiv 0$. Future SOS scripts in this repo: always
+add the $(\varepsilon, \text{cap})$ pattern. Worth a comment in
+\texttt{barrier\_sos\_2d.jl} too — it has the same vulnerability (it
+also returned $B \equiv 0$ when I re-ran it mentally).
+\end{decision}
+
+\subsection*{Open issues for follow-up}
+
+\begin{itemize}
+  \item \textbf{2-D reduction is not sound for the full plant.} The
+    dropped coupling norm is $0.459$ — non-trivial. To certify the full
+    10-state hot-standby invariance, augment the SOS state with at least
+    $T_f$ (3-D) and re-run with appropriate degree. CSDP scales
+    polynomially with monomial count; degree 4 in 3 variables is
+    $\binom{7}{3} = 35$ monomials, still very tractable.
+  \item \textbf{No disturbance.} $Q_{\mathrm{sg}} \equiv 0$ here. For
+    hot-standby with auxiliary cooling there's a residual sink. Add
+    $B_w w$ to the Lie-derivative inequality and check the worst case
+    over $w \in [-Q_{\max}, Q_{\max}]$ — standard Putinar-style
+    extension.
+  \item \textbf{Equilibrium is parametrized by IC.} Adiabatic + constant
+    rod = no unique attractor. The barrier proves invariance \emph{around}
+    $x_{\mathrm{eq}}$ but the ``correct'' $x_{\mathrm{eq}}$ depends on
+    where the system started. For a shutdown controller that actually
+    drives to a setpoint, redesign as a temperature-tracking law (PI on
+    $T_c$ with rod motion). Flagged for a separate task — current
+    \texttt{ctrl\_shutdown} is open-loop rod-held, which is honest about
+    what real plants do during a controlled hold but doesn't compose
+    with reach-avoid framing as cleanly as a tracking law would.
+  \item \textbf{Lie-derivative condition is global, not boundary.} The
+    SOS uses $-\nabla B \cdot f \in \mathrm{SOS}$ everywhere, which is
+    stronger than the Prajna $\{B = 0\}$-only condition. Convex but
+    conservative. The bilinear Putinar form
+    $-\nabla B \cdot f - \sigma_b B \in \mathrm{SOS}$ is what we'd need for
+    the boundary-only version; that's a BMI and needs alternation.
+    Acceptable conservatism here because the system has a real Lyapunov
+    function anyway.
+\end{itemize}
+
+\subsection*{Files touched}
+
+\begin{itemize}
+  \item \texttt{reachability/predicates.json} --- added \texttt{shutdown\_margin}
+    in \texttt{safety\_limits}; updated \texttt{q\_scram.X\_exit\_predicate}
+    and \texttt{X\_safe\_predicate}; left \texttt{\_X\_exit\_history} for
+    forensics.
+  \item \texttt{code/scripts/reach/reach\_scram\_pj.jl} --- added
+    \texttt{RHO\_SDM}, \texttt{SDM\_RHS}; per-horizon $\rho$-bounds + halfspace
+    LHS sup; \texttt{.mat} output extended.
+  \item \texttt{code/scripts/barrier/barrier\_sos\_2d\_shutdown.jl} ---
+    new file. Equilibrium-finder + linearization + SOS barrier.
+  \item \texttt{claude\_memory/2026-04-27-scram-X\_exit-shutdown-margin.md}
+    --- session note for the scram side.
+\end{itemize}

journal/journal.tex

@@ -76,5 +76,7 @@ Each limitation ties to a plan or an open question.
 \input{entries/2026-04-20-overnight-prompt-jump.tex}
 \newpage
 \input{entries/2026-04-21-polytopic-sos-tikhonov.tex}
+\newpage
+\input{entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex}
 
 \end{document}

presentations/prelim-presentation/outline.md

@@ -3,23 +3,23 @@
 **Format:** Assertion-evidence (Alley). Each slide: one declarative sentence
 at the top, one piece of visual evidence below. No bullet soup.
 
-**Duration:** 20 minutes. ~12 slides at ~1.5 min each (60 s typical, up to
+**Duration:** 20 minutes. ~11 content slides + title + Q&A. ~1.5–2 min/slide,
-2 min for the heavier ones).
+heavier on slides 1, 9, 10. Buffer ~1 min.
 
 **Audience:** OT-informed cybersecurity experts. Mostly CS, some control-theory
 familiarity, very little reactor-physics background. Can assume fluency with:
 LTL, automata, model-checking, reachability (as a concept), SMT/SAT. Must
 explain: PKE, reactivity, precursors, hybrid systems.
 
-**Running thread:** FRET natural-language requirement → LTL → synthesized
+**Running thread:** procedures → FRET → LTL/AIGER → DRC → continuous gotcha →
-discrete controller → continuous plant → per-mode reach-avoid proof
+plant model → reach with stiffness wall → prompt-jump fix (with validity-as-an-
-→ hybrid correctness by composition.
+invariant) → two sound nonlinear results → seam.
 
 **Design principles:**
 - **Plots over bullets.** Every result slide anchors on one figure.
-- **Physical intuition before math.** Reactor basics before PKE.
+- **Physical intuition before math.** Reactor basics in passing, not as a tutorial.
-- **Honest limitations boxed on each result slide.** Audience are cyber
+- **Honest limitations on each result slide.** Audience are cyber folks —
-  folks — they respect limits more than triumphs.
+  they respect limits more than triumphs.
 - **CS vocabulary by default, engineering terms defined inline.**
 - **End with the seam**, not with a victory lap. The thesis question is
   "how do discrete proofs and continuous proofs compose?" not "we verified
@@ -27,369 +27,462 @@ discrete controller → continuous plant → per-mode reach-avoid proof
 
 ---
 
-## Slide 1 — Title + hook
+## Slide 1 — Title + hook + scopes of control (2-frame animation)
 
-**Assertion:** Verified hybrid control for nuclear startup is a real problem
+**Assertion (frame 1):** Nuclear reactor operation is dominated by humans
-with no deployed solution.
+following procedures, and the procedures themselves have no formal verification.
 
-**Evidence:** Title block + a single image showing a control-room board of
+**Evidence (frame 1):** Full-bleed photograph of a control room — analog
-old-school analog gauges adjacent to a digital SCADA screen. Caption: "The
+gauges, paper procedures on the desk, two operators.
-gap we're filling."
 
-**Speak:** Introduce self + advisor + NRC fellowship context. Name the problem:
+**Speak (frame 1):** Self + advisor + NRC fellowship. The hook: most plants
-nuclear reactor operations are 70-80% human-error-driven at severe-accident
+are run by humans following paper procedures. We've been engineering humans
-level (Chernobyl, TMI, Fukushima all primarily human factors). Plants run on
+*out* of the loop for forty years by making the procedures more and more
-paper procedures that have no formal verification. The most advanced work to
+prescriptive — but the procedures themselves are written natural language.
-date (HARDENS, NRC-funded) verified the discrete trip system in isolation and
+We rely on humans to follow them faithfully and on tradition to keep them
-stopped there. This talk: a preliminary example of bridging discrete
+correct. (Soften: "running a nuclear reactor is well-understood" — the
-requirements all the way to continuous-state safety proofs, with the seam
+*procedures* are a knowledge-engineering artifact built over decades.)
-between them as the research contribution.
+
+**Assertion (frame 2):** Plant control decomposes into three scopes —
+strategic, operational, tactical — and formal methods most naturally land
+on the *operational* scope.
+
+**Evidence (frame 2):** Photo slides left to make room for a 3-tier pyramid
+on the right:
+
+```
+            ┌──────────────────┐
+            │   STRATEGIC      │  start it up / shut it down (hours-days)
+            └──────────────────┘
+          ┌──────────────────────┐
+          │     OPERATIONAL       │  heat up / run / scram (minutes-hours)
+          └──────────────────────┘
+        ┌──────────────────────────┐
+        │       TACTICAL            │  rod motion, valve actuation (seconds)
+        └──────────────────────────┘
+```
+
+**Speak (frame 2):** Strategic = mission-level decisions, operator authority.
+Tactical = millisecond-scale rod and valve commands, classical control. The
+**operational** middle is where the procedures live — "if cold and ready,
+heat up; if at temperature and critical, switch to power operation; if any
+trip condition, scram." This middle is where formal methods can do their
+best work because the dynamics are slow enough to verify and the logic is
+discrete enough to specify. The thesis lives here.
 
 **Reference:** thesis §2, ¶1-4. Cites Kemeny1979, Hogberg_2013, Kiniry2024.
 
-**Figures to make:** find or compose the control-room image.
+**Figures to make:** control-room photo (license-clean source TBD); pyramid
+diagram (Tikz, simple).
 
 ---
 
-## Slide 2 — The hybrid-systems gap
+## Slide 2 — FRET requirements: capture the procedure, find the seam
 
-**Assertion:** Existing tools verify either discrete *or* continuous systems,
+**Assertion:** FRET turns natural-language operational procedures into LTL,
-never the seam between them.
+and reveals the seam where a discrete predicate lands on a continuous state.
 
-**Evidence:** Two-column diagram. Left column: discrete (FRET, Spot,
+**Evidence:** A two-row figure. Top row: a FRET requirement card.
-SAL, HARDENS, TLA+). Right column: continuous (MATLAB, Modelica, CORA,
-Flow*, HyTech). A dashed line between them labeled "the bridge problem."
 
-**Speak:** HARDENS got to TRL 2-3 on the discrete side alone. Continuous-side
+> **PWR-2001:** Upon `control_mode = q_heatup ∧ t_avg_in_range ∧ p_above_crit ∧ inv1_holds`, DRC shall at the next timepoint satisfy `control_mode = q_operation`.
-reach tools exist but assume the discrete controller is given as input. What
-nobody has done: produce an *end-to-end* proof where the discrete controller's
-transitions actually fire in physical time on the modeled plant.
 
-**Reference:** thesis §2.2 "Formal Methods", §2.2.1 HARDENS, §2.2.2 Temporal
+Bottom row: the corresponding LTL, with one of the boolean atoms
-Logic. Thesis §2.3 (if exists — continuous methods).
+(`t_avg_in_range`) circled in red. Side annotation:
+
+> `t_avg_in_range` ≡ |T_c − T_c0| ≤ 2.78 °C
+> *(this is a half-space on a continuous state)*
+
+**Speak:** FRET writes requirements like "if A then next B" in a
+restricted natural-language template. Spot's compiler (FRET → LTL) checks
+that the requirement set is *realizable* — there exists a discrete
+controller that satisfies all requirements. **Here's the part that
+matters for hybrid systems:** some of these boolean atoms aren't really
+boolean. `t_avg_in_range` is a halfspace on a continuous state vector.
+The discrete controller treats it as true-or-false; the continuous
+plant has to actually *make it true* under whatever dynamics apply.
+**That gap — between the discrete requirement and the continuous
+truth-maker — is the seam.** The whole rest of the talk is about
+discharging that seam.
+
+**Reference:** `fret-pipeline/pwr_hybrid_3.json`, `reachability/predicates.json`.
+
+**Figures to make:** FRET-card visual + LTL panel.
 
 ---
 
-## Slide 3 — Running example: the PWR_HYBRID_3 controller
+## Slide 3 — From LTL to a synthesized discrete controller
 
-**Assertion:** Our target is a 4-mode hybrid controller that takes a PWR
+**Assertion:** Reactive synthesis (Spot/ltlsynt) compiles the LTL into an
-from hot standby to full power and back.
+AIGER circuit — the minimal correct discrete controller — automatically.
 
-**Evidence:** The DRC state diagram (from `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`).
+**Evidence:** Pipeline figure. Three boxes left-to-right:
-Four nodes — shutdown, heatup, operation, scram — with labeled transitions.
 
-**Speak:** One sentence of nuclear physics: a PWR has neutron population
+```
-(power), coolant temperatures at three locations, and precursor concentrations
+[ FRET requirements ]  →  [ LTL formula ]  →  [ AIGER circuit ]
-that determine delayed-neutron generation. Startup is shutdown → heatup →
+   (realizability check)    (ltlsynt synthesis)   (.aag file)
-operation; any fault triggers scram. Mode transitions gated by boolean
+```
-conditions like `t_avg_in_range ∧ p_above_crit ∧ inv1_holds`.
 
-Four modes means four *physical* behaviors the plant has to exhibit for
+Below the third box: a thumbnail of the synthesized DRC state diagram
-the discrete logic to be sound. This example stresses every layer of the
+(`fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`).
-bridge we're building.
 
-**Evidence path:** `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`,
+**Speak:** Two distinct things are happening. **FRET's realizability
-`fret-pipeline/pwr_hybrid_3.json` for one example requirement.
+check** says "your requirements are mutually consistent and a controller
+exists." **Spot/ltlsynt's reactive synthesis** actually *builds* that
+controller — it solves a parity game on the LTL formula and emits an AIGER
+circuit, the minimal-state controller satisfying the spec. We then
+extract the state diagram from the AIGER. This is well-established
+machinery in the formal-methods world; our contribution is *applying it
+to reactor operating procedures*, which is a formal-methods-free domain
+historically.
+
+**Reference:** `fret-pipeline/scripts/fret_to_synth.py`, `circuits/PWR_HYBRID_3_DRC.aag`.
 
 ---
 
-## Slide 4 — Layer 1: FRET → LTL → synthesized DRC
+## Slide 4 — The synthesized DRC for PWR_HYBRID_3
 
-**Assertion:** Natural-language requirements become a verified discrete
+**Assertion:** The discrete controller for our running example has four
-controller automatically.
+modes and seven transitions, all driven by predicates on the continuous state.
 
-**Evidence:** One-panel figure: left half = snippet of FRET natural language
+**Evidence:** Full-slide DRC state diagram from
-("Upon `control_mode = q_heatup ∧ t_avg_in_range ∧ p_above_crit ∧ inv1_holds`
+`fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`. Annotated with
-DRC shall at the next timepoint satisfy `control_mode = q_operation`"), right
+transition guards, e.g. `t_avg_in_range ∧ p_above_crit ∧ inv1_holds`
-half = the corresponding LTL, bottom = arrow labeled `ltlsynt` pointing to
+on the heatup→operation arrow.
-the synthesized state-diagram node.
 
-**Speak:** FRET compiles natural-language requirements into LTL (Spot
+**Speak:** Cold shutdown, heatup, power operation, scram. Every transition
-ecosystem). `ltlsynt` turns LTL into an AIGER circuit representing the
+is gated by a conjunction of atomic predicates. Each predicate is a
-minimal correct discrete controller. This is well-trodden ground in CS-land;
+halfspace on the 10-dimensional plant state. So the discrete controller's
-our contribution is *using it* on reactor operating procedures — so far
+correctness *as a hybrid system* depends entirely on whether the
-a formal-methods-free domain.
+continuous plant trajectory makes the right predicates true at the right
+moments. Which brings us to the gotcha.
 
-**Evidence path:** `fret-pipeline/README.md`, `pwr_hybrid_3.json` sample
+**Reference:** `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`.
-requirement, `fret-pipeline/circuits/PWR_HYBRID_3_DRC.aag` AIGER circuit.
 
 ---
 
-## Slide 5 — Layer 2: the plant model
+## Slide 5 — The continuous gotcha
 
-**Assertion:** A 10-state PKE + thermal-hydraulics model is faithful enough
+**Assertion:** A correct discrete controller does not imply a correct
-for reach analysis and tractable enough to actually compute with.
+hybrid system; the continuous-side predicates have to be discharged separately.
 
-**Evidence:** The state-vector diagram: `x = [n, C_1, ..., C_6, T_f, T_c, T_cold]`
+**Evidence:** A simple 2-panel cartoon. Left: DRC state diagram with one
-with arrows showing the coupling — neutron balance → fuel heating → coolant
+transition arrow highlighted. Right: a 2D phase portrait sketch with the
-transport → feedback to reactivity.
+trajectory drifting *outside* the predicate region — controller fires the
+transition based on logic, but the plant isn't actually where the logic
+thinks it is.
 
-**Speak:** Point-kinetic equations + three thermal nodes. 10 states, 1
+**Speak:** The DRC is correct as a discrete object — it satisfies all the
-control input (rod worth u), 1 disturbance (steam-generator heat removal Q_sg).
+LTL requirements *given* its predicate inputs. But predicates like
-Stiff system: prompt-neutron timescale is 10⁻⁴ s, thermal timescales are
+`t_avg_in_range` are continuous-state halfspaces. If the plant's
-seconds to minutes. That stiffness becomes critical later — it's what forces
+actual trajectory leaves that halfspace while the controller still
-the singular-perturbation move.
+believes it's inside, the discrete proof is meaningless. We need a
-
+*continuous-side proof* that the plant actually inhabits the right
-Five controllers: one per DRC mode (shutdown, heatup, operation under P, under
+halfspaces at the right times. That proof is per-mode and the
-LQR, scram). All Julia, all in `code/src/` + `code/controllers/`.
+methodology contribution of the chapter.
-
-**Evidence path:** `code/src/pke_th_rhs.jl`, `journal/journal.pdf` entry
-2026-04-17 for derivation.
-
-**Figures to make:** state-vector coupling diagram (can be ASCII → inkscape).
 
 ---
 
-## Slide 6 — Layer 3: the reach-avoid obligation per mode
+## Slide 6 — The plant model: 10-state PKE + thermal-hydraulics
 
-**Assertion:** Discrete correctness transfers to continuous correctness
+**Assertion:** A 10-state point-kinetic + lumped thermal-hydraulics PWR
-via one reach-avoid obligation per mode.
+model is the continuous-side surrogate — faithful enough to verify, stiff
+enough to give us trouble.
 
-**Evidence:** Taxonomy table from `reachability/WALKTHROUGH.md`:
+**Evidence:** State-vector coupling diagram showing the chain
+`u (rods) → ρ → n → P → T_f → T_c → T_cold → ρ_feedback`, with
+`Q_sg` as a bounded disturbance on `T_cold`.
+
+```
+state x = [ n  C₁..C₆  T_f  T_c  T_cold ]
+            └──prompt──┘  └─────thermal─────┘
+            (Λ ≈ 10⁻⁴ s)   (τ ≈ 10–100 s)
+                 ↑
+            stiffness ratio ≈ 10⁵
+```
+
+**Speak:** Point-kinetic equations for neutron population and six
+delayed-neutron precursor groups. Three thermal nodes — fuel, average
+coolant, cold-leg. One control input (rod-induced reactivity). One
+bounded disturbance (steam-generator heat removal). **Stiffness ratio
+of 10⁵ between prompt-neutron and thermal timescales.** Flag this now —
+it's about to be the cliff.
+
+**Reference:** `code/src/pke_th_rhs.jl`, journal entry 2026-04-17.
+
+**Figures to make:** state-vector coupling diagram (Tikz).
+
+---
+
+## Slide 7 — Per-mode reach-avoid obligations
+
+**Assertion:** Discrete-to-continuous correctness reduces to one
+reach-avoid obligation per mode — equilibrium modes prove forever-invariance,
+transition modes prove bounded-time reach-avoid.
+
+**Evidence:** Compact taxonomy:
 
 | Mode | Kind | Obligation |
 |---|---|---|
-| shutdown | equilibrium | stay in X_safe forever |
+| shutdown / operation | equilibrium | stay in safe set forever |
-| heatup | transition | from X_entry, reach X_exit in [T_min, T_max] while maintaining inv1 |
+| heatup / scram | transition | from X_entry, reach X_exit in [T_min, T_max], stay safe |
-| operation | equilibrium | stay in X_safe forever under bounded Q_sg |
-| scram | transition | from any mode, drive n safely subcritical in T_max |
 
-**Speak:** This is the structural insight. Equilibrium modes → forever-invariance
+Below the table: a tiny phase-portrait pictogram for each kind — bowl
-proofs. Transition modes → reach-avoid proofs with time budgets. The DRC
+(equilibrium) vs corridor (transition).
-transitions fire iff the reach set enters the right exit predicate. If
-every mode's obligation is discharged, the hybrid system is correct by
-composition. The methodological contribution of the chapter.
 
-**Reference:** `reachability/WALKTHROUGH.md` §1, `reachability/predicates.json`
+**Speak:** This is the structural insight. If every mode's obligation
-→ `mode_boundaries`.
+is discharged, the hybrid system is correct by composition — the discrete
+controller's transitions fire correctly because the continuous plant
+*actually arrives* at each `X_exit` predicate by the time the transition
+guard is checked. Two flavors of obligation. Two flavors of proof.
+
+**Reference:** `reachability/WALKTHROUGH.md` §1.
 
 ---
 
-## Slide 7 — Operation mode: discharge all 6 safety halfspaces (the win)
+## Slide 8 — First reach attempt: stiffness wall
 
-**Assertion:** Under LQR feedback and ±15% Q_sg variation, the operation-mode
+**Assertion:** Naive nonlinear reach (TMJets) caps out at ~10 seconds
-reach tube discharges every safety halfspace with orders of magnitude of margin.
+of horizon on the full 10-state model — orders of magnitude short of what
+the obligations need.
 
-**Evidence:** The two-panel figure `docs/figures/reach_operation_tubes.png`.
+**Evidence:** A two-panel figure. Left: a graph of horizon achieved vs
-Left: temperature tubes (T_c, T_hot, T_cold) overlaid; the near-zero width
+wall-clock time, asymptoting at ~10 s of plant time. Right: a one-line
-of T_c visually demonstrates LQR contraction. Right: ΔT_core tube with a
+caption — "T_max(heatup) = 5 hr; T_max(scram) = 60 s. We're 1800× short
-right-axis in MW showing the power is tight around nominal.
+on heatup and 6× short on scram."
 
-Per-halfspace margin table (inset or second slide if space allows):
+**Speak:** TMJets is a Taylor-model integrator — produces *sound*
+over-approximations of the nonlinear reach set. Beautiful tool. But its
+step size is bounded by the *fastest* dynamics in the system, which here
+is the prompt-neutron timescale `Λ ≈ 10⁻⁴ s`. Even on a bounded reach
+horizon, you blow your step budget propagating Taylor models that nobody
+cares about, because nothing safety-relevant happens on that timescale.
+**The stiffness ratio kills us.** We need to remove the fast modes
+without losing soundness.
 
-| Halfspace | Limit | Reach max | Margin |
+**Reference:** `code/scripts/reach/reach_heatup_nonlinear.jl`.
-|---|---:|---:|---:|
-| fuel_centerline | 1200 °C | 328.9 | **+871 K** |
-| t_avg_high_trip | 320 °C | 308.4 | **+11.6 K** |
-| n_high_trip | 1.15 | 1.012 | **+0.14** |
-
-**Speak:** LQR contracts the regulated direction (T_c) from 0.1 K to 0.03 K
-halfwidth. Uncontrolled directions (precursors) expand but don't appear in
-any safety predicate. Tightest margin is the high-flux trip at 12% headroom.
-
-**Limitations box:** linear-reach of a nonlinear plant — *approximate, not
-sound* yet for the physical plant. That's what the next slides address.
-
-**Evidence path:** `code/scripts/reach/reach_operation.jl`,
-`docs/figures/reach_operation_tubes.png`.
 
 ---
 
-## Slide 8 — Quadratic Lyapunov barrier fails structurally
+## Slide 9 — Prompt-jump reduction with validity-as-an-invariant
 
-**Assertion:** The canonical quadratic Lyapunov barrier cannot certify this
+**Assertion:** Singular-perturbation reduction eliminates the prompt
-plant — *even with LQR inside the barrier* — because of plant anisotropy.
+neutronics, gives us 30× more reach horizon — and, critically, the
+reduction's validity condition becomes part of the safety obligation
+itself, not a separate hand-wave.
 
-**Evidence:** Bar chart comparing open-loop vs LQR-closed-loop Lyapunov bound
+**Evidence:** Two stacked panels.
-on `n_high_trip` direction: OL is 27 million × nominal, CL is 1,242 × nominal.
-Title: "LQR helps 20,000× but both bounds are physically meaningless."
 
-**Speak:** Plant has prompt timescale 10⁻⁴ s vs thermal timescale ~10 s.
+Top: the algebraic substitution.
-Lyapunov matrix P inherits that 10⁴ spread. Resulting ellipsoid stretches
+$$\frac{dn}{dt} = 0 \implies n = \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho(x)}$$
-obscenely along the fast directions when projected to physical coordinates.
-No amount of controller tuning fixes this — it's set by the plant's own
-physics. The ellipsoid geometry cannot match the slab-shaped safe region.
-**This is the motivation for polynomial (SOS) or polytopic barriers.**
 
-**Reference:** `code/scripts/barrier/barrier_compare_OL_CL.jl`, journal entry
+State drops from 10 to 9. Stiffness gone. Reach steps now bounded by
-2026-04-20 (evening).
+thermal timescale, not prompt timescale.
+
+Bottom: the validity condition + Tikhonov bound.
+- Validity: `β − ρ(x) > δ > 0` over the reach set — this is a half-space
+  in the same vocabulary as every other safety predicate.
+- Error: `|x(t) − x_PJ(t)| ≤ C·Λ = O(10⁻⁴)` (Tikhonov).
+
+**Speak:** The trick that makes this thesis-shaped: we don't *assume* the
+prompt-jump reduction is valid — we **prove** it as part of the same reach
+obligation. The halfspace `prompt_critical_margin_holds` is in the per-mode
+invariant set right alongside the safety halfspaces. If reach discharges
+the invariant, it discharges *both* safety and the approximation's
+soundness. No separate validity argument. This is the
+formal-methods-shaped move I want the audience to take home.
+
+**Reference:** `code/src/pke_th_rhs_pj.jl`, journal entry 2026-04-21
+"Tikhonov bound", `reachability/predicates.json::prompt_critical_margin_*`.
 
 ---
 
-## Slide 9 — Prompt-jump reduction makes nonlinear reach tractable
+## Slide 10 — Two sound nonlinear reach proofs
 
-**Assertion:** Singular-perturbation reduction of the fast neutronics gives
+**Assertion:** With prompt-jump, both transition modes — heatup and scram —
-us 30× more reach horizon and a rigorous O(Λ) error bound via Tikhonov.
+discharge their full safety invariants over their relevant horizons.
 
-**Evidence:** Side-by-side: left panel = validate_pj_heatup.png (empirical:
+**Evidence:** Side-by-side, two panels.
-PJ and full-state trajectories overlay perfectly over 50 minutes). Right panel
-= the Tikhonov bound written out mathematically:
-$$|x(t) - x_{\mathrm{PJ}}(t)| \leq C \cdot \Lambda = \mathcal{O}(10^{-4})$$
 
-**Speak:** Setting dn/dt = 0 and solving algebraically for n gives us
+Left panel — **heatup**: tube plot from
-`n = Λ·Σλᵢ·Cᵢ / (β - ρ)`. State drops from 10 to 9, removes the Λ⁻¹
+`docs/figures/reach_heatup_pj_tubes.png`. Caption: "300 s, all 6
-stiffness. Reach tool (TMJets) can take steps on thermal timescale instead
+`inv1_holds` halfspaces discharged, T_c stable at [281.05, 291.0] °C."
-of prompt timescale.
 
-The magic: we *prove* the PJ approximation's validity condition (`β - ρ > 0`
+Right panel — **scram**: shutdown-margin discharge from
-with margin) **as part of the same reach obligation** via the
+`results/reach_scram_pj_fat.mat`. A bar chart of |ρ| vs the 1% Δk/k
-`prompt_critical_margin_heatup` halfspace in `inv1_holds`. No hand-waving —
+floor at T = 10, 30, 60 s. Caption: "Fat entry polytope (union of all
-the reach machinery proves both the safety claim and the approximation's
+mode envelopes). |ρ| ≈ 5%. **5× the requirement at 60 s.**"
-soundness together.
 
-**Evidence path:** `code/src/pke_th_rhs_pj.jl`, `code/scripts/sim/validate_pj.jl`,
+**Speak:** **Heatup**: 12,932 reach-sets, 200 s wall, tube stable —
-`journal/` entry 2026-04-21 "Tikhonov bound".
+`T_c` envelope identical at 60 s and 300 s, meaning the controller holds
+the state inside the tube indefinitely. **Scram**: from the fat entry
+polytope (any state the plant could be in across all modes plus LOCA),
+the shutdown-margin halfspace discharges at 10, 30, *and* 60 s with five
+times the NRC tech-spec margin. Both are sound for the prompt-jump-reduced
+plant; both inherit the O(Λ) Tikhonov error to the full plant. **Two
+transition modes formally verified end-to-end.** This is the headline
+result.
+
+**Limitations box:** 300 s of a 5-hour `T_max` on heatup. Step budget
+is the wall; entry refinement is the path to hours.
+
+**Evidence path:** `code/scripts/reach/reach_heatup_pj.jl`,
+`code/scripts/reach/reach_scram_pj_fat.jl`,
+`docs/figures/reach_heatup_pj_tubes.png`,
+`results/reach_scram_pj_fat.mat`.
 
 ---
 
-## Slide 10 — First sound nonlinear heatup reach (the second win)
+## Slide 11 — Composition + impact + the open seam
 
-**Assertion:** With PJ + a tightened entry box, nonlinear Taylor-model reach
+**Assertion:** Two transition modes verified, two equilibrium modes are
-discharges all 6 `inv1_holds` safety halfspaces over the first 300 s of heatup.
+the next step — the composition story holds, the open question is well-defined.
 
-**Evidence:** The four-panel `reach_heatup_pj_tubes.png`: temperature tubes
+**Evidence:** Two-column layout.
-overlaid, ΔT_core, ρ in dollars (stays between -0.25 $ and -0.05 $ —
-sub-prompt-critical), n decaying.
 
-**Speak:** TMJets Taylor-model integration. 12,932 reach-sets, 200 s wall time.
+Left column — **what's proven** (with a green check):
-Tube **stable** — T_c envelope `[281.05, 291.0]` identical at 60 s and 300 s,
-meaning the controller holds the state inside the tube indefinitely. This
-is the first sound (modulo O(Λ) PJ error) nonlinear reach-avoid artifact
-for this plant.
 
-**Limitations box:** 300 s of a 5-hour `T_max`. Step budget caps horizon; entry
+| Mode | Status |
-refinement (Blanchini-style) is the path to hours.
-
-**Evidence path:** `code/scripts/reach/reach_heatup_pj.jl configs/heatup/tight.toml`,
-`docs/figures/reach_heatup_pj_tubes.png`.
-
----
-
-## Slide 11 — Degree-4 SOS polynomial barrier, a working proof of concept
-
-**Assertion:** A degree-4 polynomial barrier certificate — the structure that
-quadratic Lyapunov cannot provide — works on the operation-mode projection.
-
-**Evidence:** Equation block showing the symbolic polynomial
-(abbreviated — the actual 13-term polynomial from `code/scripts/barrier/barrier_sos_2d.jl`).
-Side caption: "CSDP returned OPTIMAL. Solve: ~3 seconds."
-
-**Speak:** Polynomials of degree 4 can localize to slab-shaped safety regions
-in a way degree-2 (quadratic) cannot. The Prajna-Jadbabaie formulation: find
-`B(x)` such that `B ≤ 0` on entry, `B ≥ 0` on unsafe, and `∇B·f ≤ 0` on
-`{B=0}`. Sum-of-squares programming reduces this to an SDP. 2D projection
-proof of concept; scaling to full 10-state requires bigger solver (Mosek or SCS).
-
-**Reference:** `code/scripts/barrier/barrier_sos_2d.jl`, journal entry
-2026-04-21.
-
----
-
-## Slide 12 — The hybrid correctness story, assembled
-
-**Assertion:** All four pieces — FRET/synthesis, plant model, reach tubes,
-SOS/polytopic barriers — compose into a hybrid system correctness argument.
-
-**Evidence:** The composition picture. DRC state diagram at center; each
-DRC node has an arrow to a "per-mode obligation box" labeled with its proof
-artifact (tube or barrier or certificate). Arrows between nodes = transitions,
-labeled with the predicate polytope. Outer dashed box: "hybrid system
-closed-loop safety."
-
-**Speak:** What's been built, what's been proven, what's next. Transition
-correctness between modes is the next thesis-blocking piece — each mode's
-`X_exit` has to be inside the next mode's `X_entry` (with margin). That's an
-inclusion check, not a reach — but it's the final structural piece.
-
----
-
-## Slide 13 — Limitations and next steps
-
-**Assertion:** This is a preliminary example, not a deployable system.
-
-**Evidence:** A two-column list:
-
-| What's proven | What's next |
 |---|---|
-| Operation mode safe under ±15% load (approximate, via linear reach) | Nonlinear operation reach (tight SOS barrier, LOCA disturbance) |
+| Heatup | Sound nonlinear reach, 300 s, all 6 safety halfspaces |
-| Heatup PJ reach discharges all 6 safety halfspaces for 300 s | Extend to full 5-hr `T_max`; model steam-dump disturbance |
+| Scram | Sound nonlinear reach, 60 s, shutdown margin, 5× tech-spec |
-| Scram PJ n-decay monotone over 60 s | Expand `X_entry(scram)` to union of all mode tubes + LOCA |
+| FRET → AIGER | Sound discrete controller, realizability checked |
-| Degree-4 SOS barrier on 2D projection | Full 10-state SOS barrier |
+| PJ validity | Discharged inside the same reach obligation |
-| PJ error rigorously bounded by Tikhonov | Parametric α uncertainty; DNBR correlation |
 
-**Speak:** The gap from prelim example to deployable system is well-defined:
+Right column — **what's next** (amber):
-more states, tighter bounds, real tech-spec numbers, hardware-in-the-loop.
+
-None of these is a research novelty unto itself — they're engineering.
+| Open | Path |
-The research contribution is the composition framework, demonstrated end-to-end
+|---|---|
-on a nontrivial running example. Phase 2 of the thesis is filling in the
+| Operation mode forever-invariance | Polynomial barrier certificate (SOS) on PJ dynamics |
-gaps and expanding to multiple plants.
+| Hot-standby forever-invariance | Same machinery, different equilibrium |
+| Full 5-hr heatup horizon | Entry refinement (Blanchini-style) |
+| Hardware integration | Ovation DCS, scheduled |
+
+**Speak:** Where this lands. Two transition modes, formally verified
+end-to-end — discrete controller from FRET, continuous trajectory bounded
+by sound nonlinear reach, validity of the reduction proven inside the
+same obligation. **The open piece is stability proofs for the equilibrium
+modes — operation and hot-standby.** We've started on barrier certificates
+to discharge those, and the machinery works on a 2D linearization, but
+the sound treatment on the full nonlinear plant is the next thrust.
+That's the work for the next several months. **The composition framework
+holds; we're filling in cells in the matrix.**
+
+**Cyber angle close:** Formal verification of operational procedures is
+defense-in-depth for OT. Even if an attacker bypasses the comms layer
+and injects commands, a verified DRC plus a discharged reach-avoid
+envelope **constrains what the physical plant can be made to do.** That's
+an assurance axis comms-security alone can't reach.
 
 ---
 
-## Slide 14 — Questions / acknowledgements
+## Slide 12 — Q&A / acknowledgements (backup)
 
 **Assertion:** (none — backup slide)
 
-**Evidence:** Thanks to advisor, committee, NRC fellowship. Links to:
+**Evidence:** Acknowledgements. Advisor (Cole), committee, NRC fellowship.
-repo (gitea), journal PDF, thesis in progress.
+Repo (gitea), journal PDF, thesis in progress.
 
-**Speak:** Open for questions. Expect questions on:
+**Anticipated questions:**
-- "Why not Stateflow/Simulink?" → (have answer prepared; HARDENS used
+- "Why not Stateflow/Simulink?" → tool-integration story; HARDENS used
-  Cryptol for a reason — formal tool integration matters)
+  Cryptol for the same reason.
-- "What's the relationship to LTLMoP / DragonFly?" → (survey answer)
+- "How does this interact with ML components?" → out of scope; the pitch
-- "How would this interact with ML components?" → (out of scope for now,
+  is *no ML in the safety-critical loop*.
-  the whole pitch is *no ML in safety-critical loop*)
+- "What's the threat model?" → tie back to OT audience: formal methods
-- "What's your threat model for cybersecurity?" → (tie back to OT audience —
+  guarantee the controller's logic and the physical plant's behavior;
-  formal methods guarantee the controller's logic; they do not protect the
+  they do not protect comms or implementation. Defense in depth.
-  comms layer or implementation-level vulns. Mention HARDENS' focus on
+- "Why not just do nonlinear-SOS directly?" → the thrust starts there in
-  "correctness of implementation" vs our focus on "correctness of
+  the next phase; the linearized 2D version was the proof-of-concept.
-  specification")
 
 ---
 
 ## Presentation construction notes
 
-### What to build in Beamer later
+### Slide-count vs. time budget
-- Assertion-evidence template (one sentence at top, centered figure below).
-- Consistent color coding: green for "discharged" / "proved", red for
-  "limitation" / "gap", blue for "discrete-layer", amber for "continuous-layer".
-- The composition diagram on slide 12 — this will take the most Tikz work.
 
-### Figures that need to be created or dressed up for the talk
+11 content slides + title + Q&A. Average 1.6 min/slide. Allocation:
-- Slide 1: control-room visual.
-- Slide 2: discrete/continuous tools comparison.
-- Slide 3: use `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`.
-- Slide 4: FRET → LTL → AIGER panel.
-- Slide 5: 10-state coupling diagram.
-- Slide 7, 10: reuse `docs/figures/reach_*_tubes.png`.
-- Slide 8: bar chart (new, easy matplotlib).
-- Slide 9: reuse `validate_pj_heatup.png`.
-- Slide 11: latex-typeset polynomial + CSDP log snippet.
-- Slide 12: composition diagram (Tikz, will take time).
 
-### Cybersecurity angle to emphasize
+| Slide | Min | Reason |
-- The point the OT audience will care about: this kind of verification
+|---|---|---|
-  **constrains what the controller can physically do**. Even if an attacker
+| 1 (anim) | 2.5 | Hook needs riff room |
-  gets past authentication and can inject arbitrary commands, the DRC-plus-
+| 2 (FRET seam) | 2.0 | Audience unfamiliar; this is the central insight |
-  reach-certified envelope limits how bad things can get *within the
+| 3 | 1.5 | Pipeline diagram |
-  physical plant*. That's a different assurance axis than the usual
+| 4 | 1.5 | DRC walkthrough |
-  comms-security one and complements it.
+| 5 | 1.0 | Transition slide, short |
-- Formal methods as defense-in-depth: they catch bugs *before* deployment,
+| 6 | 1.5 | Plant + stiffness flag |
-  which reduces attack surface more than any runtime defense.
+| 7 | 1.0 | Taxonomy, short |
-- The PJ reduction + Tikhonov approach might be of interest for other
+| 8 | 1.5 | Wall problem |
-  safety-critical stiff systems (power grid, aerospace).
+| 9 | 2.0 | Methodology contribution; slowest |
+| 10 | 2.0 | Headline result |
+| 11 | 1.5 | Closing seam |
+| **Total** | **18.0** | + 2 min buffer |
+
+### What to build in Beamer
+
+- Assertion-evidence template — one declarative sentence at top, centered
+  figure below, optional 2-line speaker note in footer.
+- Color coding: green = sound/proven, amber = approximate/open, red = limitation,
+  blue = discrete-layer, purple = continuous-layer.
+- 2-frame animation on slide 1 (overlay specs in Beamer).
+- Inset boxes on result slides for limitations (slide 10).
+
+### Figures that need to be created or dressed up
+
+| Slide | Figure | Source / status |
+|---|---|---|
+| 1 | Control-room photo | License-clean source TBD |
+| 1 | Strategic/operational/tactical pyramid | Tikz, fresh |
+| 2 | FRET-card + LTL panel | Inkscape, fresh |
+| 3 | FRET → LTL → AIGER pipeline | Tikz, fresh |
+| 4 | DRC state diagram | `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`, dress up |
+| 5 | DRC + drifted-trajectory cartoon | Inkscape, fresh |
+| 6 | State-vector coupling diagram | Tikz, fresh |
+| 7 | Equilibrium/transition pictograms | Tikz, simple |
+| 8 | Horizon-vs-walltime stiffness graph | Matplotlib, fresh |
+| 9 | PJ algebraic substitution + Tikhonov | LaTeX math + label |
+| 10 | Heatup tubes | `docs/figures/reach_heatup_pj_tubes.png`, dress up |
+| 10 | Scram shutdown-margin bars | Matplotlib, fresh, from `reach_scram_pj_fat.mat` |
+| 11 | Two-column matrix | Beamer table |
+
+### Cybersecurity angle to emphasize (one sentence each, distributed across slides)
+
+- Slide 1: procedures are the OT-side analog of code; both deserve formal verification.
+- Slide 2: discrete-only verification (HARDENS) leaves the seam unaddressed.
+- Slide 11 close: verified physical-plant envelope is an assurance axis
+  comms security alone cannot reach.
 
 ### Things to NOT do
+
 - Don't get lost in reactor-physics detail. One sentence per physical
-  concept; get them to the CS content fast.
+  concept; keep the audience moving.
-- Don't show code unless it's a slide about *why* the code structure
+- Don't show code unless code structure is the point. Code screenshots
-  matters. Code screenshots are terrible evidence.
+  are weak evidence.
-- Don't oversell. Honest limitations at every stage builds trust with
+- Don't oversell. Honest limitations build trust with skeptical audiences.
-  a skeptical audience.
 - Don't use more than 2 bullet points on any slide. Alley rules.
+- **Don't try to fit SOS barriers as their own slide.** They live in slide 11
+  as one closing line about what's next. Cutting them out of the talk was
+  a deliberate choice — say it once, move on.
 
 ### Timing checkpoints
-- Slide 6 (mode-obligation taxonomy) by minute 8.
+
-- Slide 10 (first nonlinear reach result) by minute 13.
+- Slide 4 (DRC) by minute 6.
-- Slide 12 (composition story) by minute 17.
+- Slide 7 (taxonomy) by minute 9.
-- Slide 13 (limitations) by minute 19.
+- Slide 9 (PJ) by minute 13.
+- Slide 10 (results) by minute 16.
+- Slide 11 (close) by minute 19.
+
+### Cuts made vs. April-21 outline (for reference)
+
+- Slide 7 "Operation reach (the win)" — cut. Linear reach is a gut-check, not
+  a headline. Demoted to one bullet on slide 11.
+- Slide 8 "Quadratic Lyapunov fails" — cut. Side-quest in this arc.
+- Slide 11 "Degree-4 SOS barrier" — cut as standalone slide. Folded into
+  slide 11 "next steps" as one line.
+- Slide 6 mode-taxonomy table — folded into slide 7 informally; modes
+  appear as concrete examples (heatup, scram) when verified, not as
+  upfront enumeration.
+- Lyapunov bar-chart, per-halfspace margin table, tools-comparison
+  diagram — all cut for time.
+- Added: scopes-of-control framing on slide 1.
+- Added: continuous-gotcha transition (slide 5) explicitly calling out the seam.
+- Added: scram reach as headline result on slide 10 (was a follow-up bullet
+  on April-21 limitations slide).

reachability/predicates.json

@@ -132,6 +132,15 @@
         { "state_index": 9, "coeff": -1.0, "rhs_expr": "-275.85" }
       ],
       "_note": "This is a controller-specific specialization. If the controller or its gain change, recompute. For LQR operation and constant-u scram the margin is trivially satisfied (much larger)."
+    },
+    "shutdown_margin": {
+      "meaning": "Total reactivity is sufficiently negative — reactor safely subcritical regardless of remaining delayed-neutron flux. The actual NRC tech-spec criterion for scram success.",
+      "concretization": "rho_total <= -rho_SDM with rho_SDM = 0.01 (1% dk/k, PWR tech-spec floor). With constant u = U_SCRAM = -8*beta = -0.05202 and rho = U_SCRAM + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0), the halfspace is: alpha_f*T_f + alpha_c*T_c <= -rho_SDM - U_SCRAM + alpha_f*T_f0 + alpha_c*T_c0. Linear halfspace in (T_f, T_c).",
+      "_derivation": "alpha_f = -2.5e-5/C; alpha_c = -1e-4/C; T_f0 = 328.35 C; T_c0 = 308.35 C (from pke_params: DT_CORE = P0/(W_M*C_C) = 36.7 C, so T_HOT0 = T_COLD0 + DT_CORE = 326.7, T_c0 = (T_HOT0+T_COLD0)/2 = 308.35, T_f0 = T_c0 + P0/HA = 328.35). Concrete RHS: -0.01 + 0.05202 - 0.008209 - 0.030835 = 0.002972. At nominal: LHS = -0.039, |rho| = 5.2% — huge margin over the 1% requirement.",
+      "halfspaces": [
+        { "row": [[8, -2.5e-5], [9, -1.0e-4]], "rhs_expr": "-rho_SDM - U_SCRAM + alpha_f*T_f0 + alpha_c*T_c0" }
+      ],
+      "_note": "Replaces 'n <= 1e-4 AND T_f <= T_f0 + 50' as the scram X_exit. The power-threshold framing was structurally awkward: (1) n is nonlinear in the PJ state (n = LAMBDA * sum(LAM_i*C_i) / (BETA - rho)), (2) the conjoined T_f bound was infeasible at 60 s under decay heat (fuel doesn't cool that fast). Reactivity is the physically correct success criterion AND linear in (T_f, T_c). Verified discharged at T = 10, 30, 60 s — see results/reach_scram_pj_result.mat."
     }
   },
 
@@ -231,11 +240,12 @@
       "obligation": "from any trip-triggering state, drive reactor to safely-subcritical within T_max",
       "X_entry_description": "any state where inv1_holds or inv2_holds fails during heatup/operation",
       "X_entry_polytope": "union of (X_operation - inv2_holds-satisfying) and (X_heatup - inv1_holds-satisfying) — for demo, take x_op",
-      "X_safe_predicate": "n is monotonically non-increasing (n'(t) <= 0); T stays bounded",
+      "X_safe_predicate": "fuel_centerline AND cold_leg_subcooled (T_f, T_cold bounds maintained throughout transit; n monotonically non-increasing under U_SCRAM is implied by physics, not enforced as halfspace)",
-      "X_exit_predicate": "n <= 1e-4 AND T_f <= T_f0 + 50 C",
+      "X_exit_predicate": "shutdown_margin",
       "T_max_seconds": 60,
       "T_min_seconds": null,
-      "_T_max_rationale": "60 s. NRC requirement typically few seconds to subcritical; 60 s is generous for our lumped model with idealized rod-insertion. Real plants: rods free-fall in ~2-3 s."
+      "_T_max_rationale": "60 s. NRC requirement typically few seconds to subcritical; 60 s is generous for our lumped model with idealized rod-insertion. Real plants: rods free-fall in ~2-3 s.",
+      "_X_exit_history": "v1: 'n <= 1e-4 AND T_f <= T_f0 + 50 C'. Replaced 2026-04-27: power threshold was nonlinear in PJ state, T_f bound was infeasible at 60 s under decay heat. shutdown_margin (rho <= -0.01) is the correct NRC tech-spec criterion and a clean halfspace in (T_f, T_c)."
     }
   },