diff --git a/app/predicate_explorer.jl b/app/predicate_explorer.jl index 913ca50..36ffbe1 100644 --- a/app/predicate_explorer.jl +++ b/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 -# so the file remains a valid standalone Julia script when Pluto isn't running. +# 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). macro bind(def, element) - #= none:1 =# quote - local iv = try - Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value - catch - b -> missing - end + #! format: off + return quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value 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 @@ -60,10 +57,12 @@ md""" # ╔═╡ c995a73a-c07e-415b-abe8-57abaf9f4b46 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))), ", "))." -end + pred_path = joinpath(@__DIR__, "..", "reachability", "predicates.json") + pred_raw = JSON.parsefile(pred_path) + 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 md""" diff --git a/claude_memory/2026-04-27-scram-X_exit-shutdown-margin.md b/claude_memory/2026-04-27-scram-X_exit-shutdown-margin.md new file mode 100644 index 0000000..39a7a13 --- /dev/null +++ b/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. diff --git a/claude_memory/2026-04-28-DICE-2026-conference-intel.md b/claude_memory/2026-04-28-DICE-2026-conference-intel.md new file mode 100644 index 0000000..f3b1812 --- /dev/null +++ b/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/) diff --git a/claude_memory/2026-04-28-path1-sos-pj-sketch.md b/claude_memory/2026-04-28-path1-sos-pj-sketch.md new file mode 100644 index 0000000..4cf86d0 --- /dev/null +++ b/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. diff --git a/code/CLAUDE.md b/code/CLAUDE.md index d9b1585..995b7be 100644 --- a/code/CLAUDE.md +++ b/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 diff --git a/code/configs/heatup/tight.toml b/code/configs/heatup/tight.toml index df3dce9..53d1401 100644 --- a/code/configs/heatup/tight.toml +++ b/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 diff --git a/code/scripts/barrier/barrier_sos_2d.jl b/code/scripts/barrier/barrier_sos_2d.jl index 7ad302c..b97da3f 100644 --- a/code/scripts/barrier/barrier_sos_2d.jl +++ b/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 -# closed-loop. If this works, scaling to full 10-state is a matter -# of increasing degree and throughput. +# barrier certificate on a reduced model. If this works, scaling to +# full 10-state is a matter of increasing degree and throughput. # -# Reduced dynamics: project the LQR closed-loop onto (dT_c, dn), the -# primary safety direction and the dominant unregulated direction. -# A_red, B_w_red are the 2x2 / 2x1 submatrices corresponding to these +# Reduced dynamics: project the LQR closed-loop onto (dn, dT_c), the +# dominant unregulated direction and the primary safety direction. +# 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. -# -# 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. +# Unsafe focus: dn ≥ +0.15 (high-flux trip; the harder direction +# under LQR because positive-n excursions trip n_high before T_c trips). 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) -println() +show(stdout, "text/plain", A_cl_red); 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 --- -# dx = [dn; dTc] = [x[1]; x[2]] in polynomial variables. -@polyvar x1 x2 +# --- SOS sets --- +@polyvar x1 x2 # x1 = dn, x2 = dT_c -# Dynamics with worst-case constant w: -w_bar = 0.15 * plant.P0 -# Split disturbance into its mid + extreme, handle as bounded constant. -# For the Lie derivative check we use the WORST-CASE w that maximizes -# the outward velocity. Since B_w_red is a known 2-vector and ∂B/∂x -# is polynomial in x, the max-over-w is achieved at w ∈ {-w_bar, +w_bar}. -# Defer that max — check both worst cases separately. +entry_halfspaces = [ + 0.01 - x1, # dn ≤ 0.01 + x1 + 0.01, # dn ≥ -0.01 + 0.1 - x2, # dT_c ≤ 0.1 + x2 + 0.1, # dT_c ≥ -0.1 +] -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: -# g1 = 5 - x2 (dT_c ≤ 5) -# g2 = x2 + 5 (dT_c ≥ -5) -# g3 = 0.15 - x1 (dn ≤ 0.15) -# g4 = x1 + 0.15 (dn ≥ -0.15) -# Unsafe set = complement; for SOS we use the Putinar formulation where -# 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. +# --- Solve --- +println(" Solving SOS feasibility (degree-4 B, ε-slack capped at 1.0)...") +result = solve_sos_barrier_2d(A_cl_red, (x1, x2), + entry_halfspaces, unsafe_halfspaces; + barrier_degree=4, multiplier_degree=2, + eps_cap=1.0) -# Entry set: -# g_e1 = 0.1 - x2; g_e2 = x2 + 0.1; g_e3 = 0.01 - x1; g_e4 = x1 + 0.01. - -g_s1 = 5.0 - x2 -g_s2 = x2 + 5.0 -g_s3 = 0.15 - x1 -g_s4 = x1 + 0.15 -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.") +println(" Status: $(result.status)") +if result.status == MOI.OPTIMAL && result.ε > 1e-8 + println(" ✅ ε* = $(round(result.ε; digits=4)) — real certificate.") + println(" B(x) = $(result.B)") +elseif result.status == MOI.OPTIMAL + println(" ⚠ ε ≈ 0 — solver returned trivial B ≡ 0. No real barrier") + println(" at degree 4 with these sets.") else - println(" ⚠ Solver stopped with: $status") + println(" ❌ $(result.status). Try higher degree or relax sets.") end diff --git a/code/scripts/barrier/barrier_sos_2d_shutdown.jl b/code/scripts/barrier/barrier_sos_2d_shutdown.jl new file mode 100644 index 0000000..efb59de --- /dev/null +++ b/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)) diff --git a/code/scripts/plot/plot_reach_tubes.jl b/code/scripts/plot/plot_reach_tubes.jl index 998107b..011b44e 100644 --- a/code/scripts/plot/plot_reach_tubes.jl +++ b/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() - mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability", - "reach_heatup_pj_tight_full.mat") +function plot_tubes_heatup_pj(probe_seconds::Int=8000) + # Reads the per-probe save format from reach_heatup_pj.jl with the + # 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"]) - Tc_lo = vec(d["Tc_lo_ts"]); Tc_hi = vec(d["Tc_hi_ts"]) - Tf_lo = vec(d["Tf_lo_ts"]); Tf_hi = vec(d["Tf_hi_ts"]) - Tco_lo = vec(d["Tco_lo_ts"]); Tco_hi = vec(d["Tco_hi_ts"]) - n_lo = vec(d["n_lo_ts"]); n_hi = vec(d["n_hi_ts"]) - rho_lo = vec(d["rho_lo_ts"]); rho_hi = vec(d["rho_hi_ts"]) + pre = "T_$(probe_seconds)_" + if !haskey(d, pre * "t_arr") + # Fall back to legacy bare-key format (older mat files). + pre = "" + haskey(d, "t_arr") || error("Neither '$pre' nor bare keys found in $mat_path. Available probes: " * + 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") - mat_path = joinpath(@__DIR__, "..", "..", "..", "reachability", - "reach_heatup_pj_tight_full.mat") +if which_plot in ("scram_full", "both") + mat_path = joinpath(@__DIR__, "..", "..", "..", "results", + "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("(Run scripts/reach_heatup_pj_tight_full.jl first.)") + println("Skipping scram_full plot — $mat_path not found.") + 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 diff --git a/code/scripts/reach/reach_heatup_pj.jl b/code/scripts/reach/reach_heatup_pj.jl index f5645f6..c040b98 100644 --- a/code/scripts/reach/reach_heatup_pj.jl +++ b/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_hi = T_STANDBY + RAMP_RATE_CS * t_hi_here + Tref_lo = T_REF_START_C + RAMP_RATE_CS * t_lo_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 diff --git a/code/scripts/reach/reach_scram_full_fat.jl b/code/scripts/reach/reach_scram_full_fat.jl new file mode 100644 index 0000000..f690bb5 --- /dev/null +++ b/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") diff --git a/code/scripts/reach/reach_scram_pj.jl b/code/scripts/reach/reach_scram_pj.jl index 9b3d4f4..ea07c10 100644 --- a/code/scripts/reach/reach_scram_pj.jl +++ b/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 -# n <= 1e-4 within T_max = 60 s. Constant control u = -8*beta (rods -# slammed in). Q_sg = 3% P0 (decay-heat-level sink, placeholder). +# Scram obligation: from any operating-envelope state, drive total +# reactivity below the shutdown-margin threshold (rho <= -0.01, i.e. +# 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" - 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") + 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 — 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") diff --git a/code/scripts/reach/reach_scram_pj_fat.jl b/code/scripts/reach/reach_scram_pj_fat.jl index ff12ade..fec3d16 100644 --- a/code/scripts/reach/reach_scram_pj_fat.jl +++ b/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) diff --git a/code/scripts/sim/main_mode_sweep.jl b/code/scripts/sim/main_mode_sweep.jl index 9a7c228..cbb924f 100644 --- a/code/scripts/sim/main_mode_sweep.jl +++ b/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) -T_standby = plant.T_c0 + pred_raw["derived"]["T_standby_offset_C"] - -# --- ICs --- -x_op = pke_initial_conditions(plant) - -# 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] +predicates = JSON.parsefile(pred_path) +states = pke_states(plant; predicates=predicates) +T_standby = states.T_standby +x_op = states.op +x_shut = states.shutdown +x_heat = states.heatup # --- LQR gain (needed by ctrl_operation_lqr_factory) --- A, B, B_w, _, _, _ = pke_linearize(plant) diff --git a/code/scripts/sim/validate_pj.jl b/code/scripts/sim/validate_pj.jl index 719bbf9..2212ab2 100644 --- a/code/scripts/sim/validate_pj.jl +++ b/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. -n0 = 1e-3 -C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n0 -x0_full = [n0; C0; T_standby; T_standby; T_standby] -x0_pj = [C0; T_standby; T_standby; T_standby] +# Initial conditions — both at n=1e-3 (the heatup canonical IC), same temperatures. +x0_full = states.heatup # 10-state +x0_pj = states.heatup[2:end] # drop n; PJ reconstructs it tspan = (0.0, 3000.0) diff --git a/code/src/pke_states.jl b/code/src/pke_states.jl new file mode 100644 index 0000000..e49d2cc --- /dev/null +++ b/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..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 diff --git a/code/src/sos_barrier.jl b/code/src/sos_barrier.jl new file mode 100644 index 0000000..dcd860d --- /dev/null +++ b/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 diff --git a/docs/figures/reach_heatup_pj_tubes.png b/docs/figures/reach_heatup_pj_tubes.png index 11aad3b..55c77f5 100644 Binary files a/docs/figures/reach_heatup_pj_tubes.png and b/docs/figures/reach_heatup_pj_tubes.png differ diff --git a/docs/figures/reach_scram_full_tubes.png b/docs/figures/reach_scram_full_tubes.png new file mode 100644 index 0000000..c60d859 Binary files /dev/null and b/docs/figures/reach_scram_full_tubes.png differ diff --git a/docs/figures/reach_scram_pj_fat_tubes.png b/docs/figures/reach_scram_pj_fat_tubes.png new file mode 100644 index 0000000..d8c5949 Binary files /dev/null and b/docs/figures/reach_scram_pj_fat_tubes.png differ diff --git a/docs/figures/reach_scram_pj_tubes.png b/docs/figures/reach_scram_pj_tubes.png new file mode 100644 index 0000000..1e4b090 Binary files /dev/null and b/docs/figures/reach_scram_pj_tubes.png differ diff --git a/docs/model_cheatsheet.md b/docs/model_cheatsheet.md new file mode 100644 index 0000000..b36266b --- /dev/null +++ b/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` | diff --git a/journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex b/journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex new file mode 100644 index 0000000..5631f03 --- /dev/null +++ b/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} diff --git a/journal/journal.tex b/journal/journal.tex index e92c53a..339c88d 100644 --- a/journal/journal.tex +++ b/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} diff --git a/presentations/prelim-presentation/outline.md b/presentations/prelim-presentation/outline.md index 001f140..95ff57b 100644 --- a/presentations/prelim-presentation/outline.md +++ b/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 -2 min for the heavier ones). +**Duration:** 20 minutes. ~11 content slides + title + Q&A. ~1.5–2 min/slide, +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 -discrete controller → continuous plant → per-mode reach-avoid proof -→ hybrid correctness by composition. +**Running thread:** procedures → FRET → LTL/AIGER → DRC → continuous gotcha → +plant model → reach with stiffness wall → prompt-jump fix (with validity-as-an- +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. -- **Honest limitations boxed on each result slide.** Audience are cyber - folks — they respect limits more than triumphs. +- **Physical intuition before math.** Reactor basics in passing, not as a tutorial. +- **Honest limitations on each result slide.** Audience are cyber folks — + 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 -with no deployed solution. +**Assertion (frame 1):** Nuclear reactor operation is dominated by humans +following procedures, and the procedures themselves have no formal verification. -**Evidence:** Title block + a single image showing a control-room board of -old-school analog gauges adjacent to a digital SCADA screen. Caption: "The -gap we're filling." +**Evidence (frame 1):** Full-bleed photograph of a control room — analog +gauges, paper procedures on the desk, two operators. -**Speak:** Introduce self + advisor + NRC fellowship context. Name the problem: -nuclear reactor operations are 70-80% human-error-driven at severe-accident -level (Chernobyl, TMI, Fukushima all primarily human factors). Plants run on -paper procedures that have no formal verification. The most advanced work to -date (HARDENS, NRC-funded) verified the discrete trip system in isolation and -stopped there. This talk: a preliminary example of bridging discrete -requirements all the way to continuous-state safety proofs, with the seam -between them as the research contribution. +**Speak (frame 1):** Self + advisor + NRC fellowship. The hook: most plants +are run by humans following paper procedures. We've been engineering humans +*out* of the loop for forty years by making the procedures more and more +prescriptive — but the procedures themselves are written natural language. +We rely on humans to follow them faithfully and on tradition to keep them +correct. (Soften: "running a nuclear reactor is well-understood" — the +*procedures* are a knowledge-engineering artifact built over decades.) + +**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, -never the seam between them. +**Assertion:** FRET turns natural-language operational procedures into LTL, +and reveals the seam where a discrete predicate lands on a continuous state. -**Evidence:** Two-column diagram. Left column: discrete (FRET, Spot, -SAL, HARDENS, TLA+). Right column: continuous (MATLAB, Modelica, CORA, -Flow*, HyTech). A dashed line between them labeled "the bridge problem." +**Evidence:** A two-row figure. Top row: a FRET requirement card. -**Speak:** HARDENS got to TRL 2-3 on the discrete side alone. Continuous-side -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. +> **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`. -**Reference:** thesis §2.2 "Formal Methods", §2.2.1 HARDENS, §2.2.2 Temporal -Logic. Thesis §2.3 (if exists — continuous methods). +Bottom row: the corresponding LTL, with one of the boolean atoms +(`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 -from hot standby to full power and back. +**Assertion:** Reactive synthesis (Spot/ltlsynt) compiles the LTL into an +AIGER circuit — the minimal correct discrete controller — automatically. -**Evidence:** The DRC state diagram (from `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`). -Four nodes — shutdown, heatup, operation, scram — with labeled transitions. +**Evidence:** Pipeline figure. Three boxes left-to-right: -**Speak:** One sentence of nuclear physics: a PWR has neutron population -(power), coolant temperatures at three locations, and precursor concentrations -that determine delayed-neutron generation. Startup is shutdown → heatup → -operation; any fault triggers scram. Mode transitions gated by boolean -conditions like `t_avg_in_range ∧ p_above_crit ∧ inv1_holds`. +``` +[ FRET requirements ] → [ LTL formula ] → [ AIGER circuit ] + (realizability check) (ltlsynt synthesis) (.aag file) +``` -Four modes means four *physical* behaviors the plant has to exhibit for -the discrete logic to be sound. This example stresses every layer of the -bridge we're building. +Below the third box: a thumbnail of the synthesized DRC state diagram +(`fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`). -**Evidence path:** `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`, -`fret-pipeline/pwr_hybrid_3.json` for one example requirement. +**Speak:** Two distinct things are happening. **FRET's realizability +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 -controller automatically. +**Assertion:** The discrete controller for our running example has four +modes and seven transitions, all driven by predicates on the continuous state. -**Evidence:** One-panel figure: left half = snippet of FRET natural language -("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`"), right -half = the corresponding LTL, bottom = arrow labeled `ltlsynt` pointing to -the synthesized state-diagram node. +**Evidence:** Full-slide DRC state diagram from +`fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`. Annotated with +transition guards, e.g. `t_avg_in_range ∧ p_above_crit ∧ inv1_holds` +on the heatup→operation arrow. -**Speak:** FRET compiles natural-language requirements into LTL (Spot -ecosystem). `ltlsynt` turns LTL into an AIGER circuit representing the -minimal correct discrete controller. This is well-trodden ground in CS-land; -our contribution is *using it* on reactor operating procedures — so far -a formal-methods-free domain. +**Speak:** Cold shutdown, heatup, power operation, scram. Every transition +is gated by a conjunction of atomic predicates. Each predicate is a +halfspace on the 10-dimensional plant state. So the discrete controller's +correctness *as a hybrid system* depends entirely on whether the +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 -requirement, `fret-pipeline/circuits/PWR_HYBRID_3_DRC.aag` AIGER circuit. +**Reference:** `fret-pipeline/diagrams/PWR_HYBRID_3_DRC_states.png`. --- -## Slide 5 — Layer 2: the plant model +## Slide 5 — The continuous gotcha -**Assertion:** A 10-state PKE + thermal-hydraulics model is faithful enough -for reach analysis and tractable enough to actually compute with. +**Assertion:** A correct discrete controller does not imply a correct +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]` -with arrows showing the coupling — neutron balance → fuel heating → coolant -transport → feedback to reactivity. +**Evidence:** A simple 2-panel cartoon. Left: DRC state diagram with one +transition arrow highlighted. Right: a 2D phase portrait sketch with the +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 -control input (rod worth u), 1 disturbance (steam-generator heat removal Q_sg). -Stiff system: prompt-neutron timescale is 10⁻⁴ s, thermal timescales are -seconds to minutes. That stiffness becomes critical later — it's what forces -the singular-perturbation move. - -Five controllers: one per DRC mode (shutdown, heatup, operation under P, under -LQR, scram). All Julia, all in `code/src/` + `code/controllers/`. - -**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). +**Speak:** The DRC is correct as a discrete object — it satisfies all the +LTL requirements *given* its predicate inputs. But predicates like +`t_avg_in_range` are continuous-state halfspaces. If the plant's +actual trajectory leaves that halfspace while the controller still +believes it's inside, the discrete proof is meaningless. We need a +*continuous-side proof* that the plant actually inhabits the right +halfspaces at the right times. That proof is per-mode and the +methodology contribution of the chapter. --- -## 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 -via one reach-avoid obligation per mode. +**Assertion:** A 10-state point-kinetic + lumped thermal-hydraulics PWR +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 | -| heatup | transition | from X_entry, reach X_exit in [T_min, T_max] while maintaining inv1 | -| operation | equilibrium | stay in X_safe forever under bounded Q_sg | -| scram | transition | from any mode, drive n safely subcritical in T_max | +| shutdown / operation | equilibrium | stay in safe set forever | +| heatup / scram | transition | from X_entry, reach X_exit in [T_min, T_max], stay safe | -**Speak:** This is the structural insight. Equilibrium modes → forever-invariance -proofs. Transition modes → reach-avoid proofs with time budgets. The DRC -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. +Below the table: a tiny phase-portrait pictogram for each kind — bowl +(equilibrium) vs corridor (transition). -**Reference:** `reachability/WALKTHROUGH.md` §1, `reachability/predicates.json` -→ `mode_boundaries`. +**Speak:** This is the structural insight. If every mode's obligation +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 -reach tube discharges every safety halfspace with orders of magnitude of margin. +**Assertion:** Naive nonlinear reach (TMJets) caps out at ~10 seconds +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`. -Left: temperature tubes (T_c, T_hot, T_cold) overlaid; the near-zero width -of T_c visually demonstrates LQR contraction. Right: ΔT_core tube with a -right-axis in MW showing the power is tight around nominal. +**Evidence:** A two-panel figure. Left: a graph of horizon achieved vs +wall-clock time, asymptoting at ~10 s of plant time. Right: a one-line +caption — "T_max(heatup) = 5 hr; T_max(scram) = 60 s. We're 1800× short +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 | -|---|---:|---:|---:| -| 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`. +**Reference:** `code/scripts/reach/reach_heatup_nonlinear.jl`. --- -## 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 -plant — *even with LQR inside the barrier* — because of plant anisotropy. +**Assertion:** Singular-perturbation reduction eliminates the prompt +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 -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." +**Evidence:** Two stacked panels. -**Speak:** Plant has prompt timescale 10⁻⁴ s vs thermal timescale ~10 s. -Lyapunov matrix P inherits that 10⁴ spread. Resulting ellipsoid stretches -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.** +Top: the algebraic substitution. +$$\frac{dn}{dt} = 0 \implies n = \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho(x)}$$ -**Reference:** `code/scripts/barrier/barrier_compare_OL_CL.jl`, journal entry -2026-04-20 (evening). +State drops from 10 to 9. Stiffness gone. Reach steps now bounded by +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 -us 30× more reach horizon and a rigorous O(Λ) error bound via Tikhonov. +**Assertion:** With prompt-jump, both transition modes — heatup and scram — +discharge their full safety invariants over their relevant horizons. -**Evidence:** Side-by-side: left panel = validate_pj_heatup.png (empirical: -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})$$ +**Evidence:** Side-by-side, two panels. -**Speak:** Setting dn/dt = 0 and solving algebraically for n gives us -`n = Λ·Σλᵢ·Cᵢ / (β - ρ)`. State drops from 10 to 9, removes the Λ⁻¹ -stiffness. Reach tool (TMJets) can take steps on thermal timescale instead -of prompt timescale. +Left panel — **heatup**: tube plot from +`docs/figures/reach_heatup_pj_tubes.png`. Caption: "300 s, all 6 +`inv1_holds` halfspaces discharged, T_c stable at [281.05, 291.0] °C." -The magic: we *prove* the PJ approximation's validity condition (`β - ρ > 0` -with margin) **as part of the same reach obligation** via the -`prompt_critical_margin_heatup` halfspace in `inv1_holds`. No hand-waving — -the reach machinery proves both the safety claim and the approximation's -soundness together. +Right panel — **scram**: shutdown-margin discharge from +`results/reach_scram_pj_fat.mat`. A bar chart of |ρ| vs the 1% Δk/k +floor at T = 10, 30, 60 s. Caption: "Fat entry polytope (union of all +mode envelopes). |ρ| ≈ 5%. **5× the requirement at 60 s.**" -**Evidence path:** `code/src/pke_th_rhs_pj.jl`, `code/scripts/sim/validate_pj.jl`, -`journal/` entry 2026-04-21 "Tikhonov bound". +**Speak:** **Heatup**: 12,932 reach-sets, 200 s wall, tube stable — +`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 -discharges all 6 `inv1_holds` safety halfspaces over the first 300 s of heatup. +**Assertion:** Two transition modes verified, two equilibrium modes are +the next step — the composition story holds, the open question is well-defined. -**Evidence:** The four-panel `reach_heatup_pj_tubes.png`: temperature tubes -overlaid, ΔT_core, ρ in dollars (stays between -0.25 $ and -0.05 $ — -sub-prompt-critical), n decaying. +**Evidence:** Two-column layout. -**Speak:** TMJets Taylor-model integration. 12,932 reach-sets, 200 s wall time. -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. +Left column — **what's proven** (with a green check): -**Limitations box:** 300 s of a 5-hour `T_max`. Step budget caps horizon; entry -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 | +| Mode | Status | |---|---| -| Operation mode safe under ±15% load (approximate, via linear reach) | Nonlinear operation reach (tight SOS barrier, LOCA disturbance) | -| Heatup PJ reach discharges all 6 safety halfspaces for 300 s | Extend to full 5-hr `T_max`; model steam-dump disturbance | -| Scram PJ n-decay monotone over 60 s | Expand `X_entry(scram)` to union of all mode tubes + LOCA | -| Degree-4 SOS barrier on 2D projection | Full 10-state SOS barrier | -| PJ error rigorously bounded by Tikhonov | Parametric α uncertainty; DNBR correlation | +| Heatup | Sound nonlinear reach, 300 s, all 6 safety halfspaces | +| Scram | Sound nonlinear reach, 60 s, shutdown margin, 5× tech-spec | +| FRET → AIGER | Sound discrete controller, realizability checked | +| PJ validity | Discharged inside the same reach obligation | -**Speak:** The gap from prelim example to deployable system is well-defined: -more states, tighter bounds, real tech-spec numbers, hardware-in-the-loop. -None of these is a research novelty unto itself — they're engineering. -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 -gaps and expanding to multiple plants. +Right column — **what's next** (amber): + +| Open | Path | +|---|---| +| Operation mode forever-invariance | Polynomial barrier certificate (SOS) on PJ dynamics | +| 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: -repo (gitea), journal PDF, thesis in progress. +**Evidence:** Acknowledgements. Advisor (Cole), committee, NRC fellowship. +Repo (gitea), journal PDF, thesis in progress. -**Speak:** Open for questions. Expect questions on: -- "Why not Stateflow/Simulink?" → (have answer prepared; HARDENS used - Cryptol for a reason — formal tool integration matters) -- "What's the relationship to LTLMoP / DragonFly?" → (survey answer) -- "How would this interact with ML components?" → (out of scope for now, - the whole pitch is *no ML in safety-critical loop*) -- "What's your threat model for cybersecurity?" → (tie back to OT audience — - formal methods guarantee the controller's logic; they do not protect the - comms layer or implementation-level vulns. Mention HARDENS' focus on - "correctness of implementation" vs our focus on "correctness of - specification") +**Anticipated questions:** +- "Why not Stateflow/Simulink?" → tool-integration story; HARDENS used + Cryptol for the same reason. +- "How does this interact with ML components?" → out of scope; the pitch + is *no ML in the safety-critical loop*. +- "What's the threat model?" → tie back to OT audience: formal methods + guarantee the controller's logic and the physical plant's behavior; + they do not protect comms or implementation. Defense in depth. +- "Why not just do nonlinear-SOS directly?" → the thrust starts there in + the next phase; the linearized 2D version was the proof-of-concept. --- ## Presentation construction notes -### What to build in Beamer later -- 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. +### Slide-count vs. time budget -### Figures that need to be created or dressed up for the talk -- 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). +11 content slides + title + Q&A. Average 1.6 min/slide. Allocation: -### Cybersecurity angle to emphasize -- The point the OT audience will care about: this kind of verification - **constrains what the controller can physically do**. Even if an attacker - gets past authentication and can inject arbitrary commands, the DRC-plus- - reach-certified envelope limits how bad things can get *within the - physical plant*. That's a different assurance axis than the usual - comms-security one and complements it. -- Formal methods as defense-in-depth: they catch bugs *before* deployment, - which reduces attack surface more than any runtime defense. -- The PJ reduction + Tikhonov approach might be of interest for other - safety-critical stiff systems (power grid, aerospace). +| Slide | Min | Reason | +|---|---|---| +| 1 (anim) | 2.5 | Hook needs riff room | +| 2 (FRET seam) | 2.0 | Audience unfamiliar; this is the central insight | +| 3 | 1.5 | Pipeline diagram | +| 4 | 1.5 | DRC walkthrough | +| 5 | 1.0 | Transition slide, short | +| 6 | 1.5 | Plant + stiffness flag | +| 7 | 1.0 | Taxonomy, short | +| 8 | 1.5 | Wall problem | +| 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. -- Don't show code unless it's a slide about *why* the code structure - matters. Code screenshots are terrible evidence. -- Don't oversell. Honest limitations at every stage builds trust with - a skeptical audience. + concept; keep the audience moving. +- Don't show code unless code structure is the point. Code screenshots + are weak evidence. +- Don't oversell. Honest limitations build trust with skeptical audiences. - 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 12 (composition story) by minute 17. -- Slide 13 (limitations) by minute 19. + +- Slide 4 (DRC) by minute 6. +- Slide 7 (taxonomy) by minute 9. +- 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). diff --git a/reachability/predicates.json b/reachability/predicates.json index 3dff282..bd19716 100644 --- a/reachability/predicates.json +++ b/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_exit_predicate": "n <= 1e-4 AND T_f <= T_f0 + 50 C", + "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": "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)." } },