Compare commits
1 Commits
main
...
draft/2026
| Author | SHA1 | Date | |
|---|---|---|---|
|
c5133401e0 |
@ -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"""
|
||||
|
||||
88
claude_memory/2026-04-27-scram-X_exit-shutdown-margin.md
Normal file
88
claude_memory/2026-04-27-scram-X_exit-shutdown-margin.md
Normal file
@ -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.
|
||||
268
claude_memory/2026-04-28-DICE-2026-conference-intel.md
Normal file
268
claude_memory/2026-04-28-DICE-2026-conference-intel.md
Normal file
@ -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/)
|
||||
125
claude_memory/2026-04-28-path1-sos-pj-sketch.md
Normal file
125
claude_memory/2026-04-28-path1-sos-pj-sketch.md
Normal file
@ -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.
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
119
code/scripts/barrier/barrier_sos_2d_shutdown.jl
Normal file
119
code/scripts/barrier/barrier_sos_2d_shutdown.jl
Normal file
@ -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))
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
284
code/scripts/reach/reach_scram_full_fat.jl
Normal file
284
code/scripts/reach/reach_scram_full_fat.jl
Normal file
@ -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")
|
||||
@ -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")
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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)
|
||||
|
||||
|
||||
61
code/src/pke_states.jl
Normal file
61
code/src/pke_states.jl
Normal file
@ -0,0 +1,61 @@
|
||||
"""
|
||||
pke_states(plant; predicates=nothing,
|
||||
T_standby_offset_C=-33.333333,
|
||||
n_shutdown=1e-6, n_heatup=1e-3)
|
||||
|
||||
Single-point initial conditions for every canonical DRC mode.
|
||||
|
||||
Returns a NamedTuple with:
|
||||
|
||||
- `T_standby` — scalar, hot-standby coolant temperature [°C]
|
||||
- `op` — full-power steady state (`n = 1`, precursors at equilibrium)
|
||||
- `shutdown` — hot standby (trace `n`, all temperatures at `T_standby`)
|
||||
- `heatup` — critical at low power, thermals matched to `T_standby`
|
||||
|
||||
Each state is a 10-element `Vector{Float64}` matching
|
||||
`[n; C₁..C₆; T_f; T_c; T_cold]`.
|
||||
|
||||
# Source of truth for `T_standby`
|
||||
|
||||
Pass `predicates=JSON.parsefile("reachability/predicates.json")` to read
|
||||
the canonical offset from `predicates.derived.T_standby_offset_C`. If
|
||||
omitted, falls back to the keyword default (`-33.333333` °C, equivalent
|
||||
to `-60` °F — matches the `predicates.json` value at the time of writing).
|
||||
|
||||
# Why a single function
|
||||
|
||||
ICs were previously duplicated across `main_mode_sweep.jl`,
|
||||
`barrier_sos_2d_shutdown.jl`, `validate_pj.jl`, etc. Each copy independently
|
||||
defined `T_standby`, `n_shut = 1e-6`, `n_heat = 1e-3`, and the precursor-
|
||||
equilibrium formula. Lifting once means you change power-level conventions
|
||||
in one place.
|
||||
|
||||
# X_entry polytopes
|
||||
|
||||
This function returns single-point ICs for forward simulation and
|
||||
equilibrium finding. The `X_entry` *polytopes* used by reach analysis
|
||||
live in `reachability/predicates.json::mode_definitions.<mode>.X_entry_polytope`
|
||||
and should be loaded directly by reach scripts — don't conflate the two.
|
||||
"""
|
||||
function pke_states(plant; predicates=nothing,
|
||||
T_standby_offset_C=-33.333333,
|
||||
n_shutdown=1e-6, n_heatup=1e-3)
|
||||
if predicates !== nothing
|
||||
T_standby = plant.T_c0 + predicates["derived"]["T_standby_offset_C"]
|
||||
else
|
||||
T_standby = plant.T_c0 + T_standby_offset_C
|
||||
end
|
||||
|
||||
# Full-power steady state.
|
||||
op = pke_initial_conditions(plant)
|
||||
|
||||
# Hot-standby: trace n, all thermals at T_standby.
|
||||
C_shut = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_shutdown
|
||||
shutdown = [n_shutdown; C_shut; T_standby; T_standby; T_standby]
|
||||
|
||||
# Heatup: critical at 0.1% power, thermals matched to T_standby.
|
||||
C_heat = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_heatup
|
||||
heatup = [n_heatup; C_heat; T_standby; T_standby; T_standby]
|
||||
|
||||
return (; T_standby, op, shutdown, heatup)
|
||||
end
|
||||
118
code/src/sos_barrier.jl
Normal file
118
code/src/sos_barrier.jl
Normal file
@ -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
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 90 KiB After Width: | Height: | Size: 109 KiB |
BIN
docs/figures/reach_scram_full_tubes.png
Normal file
BIN
docs/figures/reach_scram_full_tubes.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 136 KiB |
BIN
docs/figures/reach_scram_pj_fat_tubes.png
Normal file
BIN
docs/figures/reach_scram_pj_fat_tubes.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 118 KiB |
BIN
docs/figures/reach_scram_pj_tubes.png
Normal file
BIN
docs/figures/reach_scram_pj_tubes.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 115 KiB |
220
docs/model_cheatsheet.md
Normal file
220
docs/model_cheatsheet.md
Normal file
@ -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` |
|
||||
217
journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex
Normal file
217
journal/entries/2026-04-27-shutdown-sos-and-scram-X_exit.tex
Normal file
@ -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}
|
||||
@ -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}
|
||||
|
||||
@ -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).
|
||||
|
||||
@ -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)."
|
||||
}
|
||||
},
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user