diff --git a/CLAUDE.md b/CLAUDE.md index 68090fc..2c9ca8a 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -41,8 +41,9 @@ hybrid controller. This repo composes three layers of the same story: |---|---|---|---| | Story | `thesis/` (submodule) | The PhD proposal: math, prose, claims | Motivates and defines the methodology | | Discrete | `fret-pipeline/` | FRET requirements → LTL → ltlsynt → AIGER circuit → state machine | Implements Thrust 1 + 2 | -| Continuous | `plant-model/` | 10-state PKE + thermal-hydraulics PWR model (MATLAB/Octave) | Plant for Thrust 3 (reachability / barrier certs) | -| Verification | `reachability/` | Continuous-mode verification artifacts | Thrust 3 output — currently empty | +| Continuous + Verification | `code/` + `reachability/` | 10-state PKE + thermal-hydraulics PWR model, controllers, reach machinery (all Julia) plus `predicates.json` and `WALKTHROUGH.md` | Plant + Thrust 3 (reachability / barrier certs) | +| Journal | `journal/` | LaTeX invention log, dated entries | Project history / re-derivability | +| App | `app/` | Pluto.jl predicate explorer | FRET-adjacent UI for hybrid-systems groups | | Hardware | `hardware/` | Ovation HIL artifacts | Integration milestone — currently empty | | Deliverables | `presentations/` | Conference papers and talks (submodules) | Dissemination | @@ -57,24 +58,31 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM). `synthesize.sh` → re-trace via `trace_aiger.py`. New `.svg` / `.dot` / `.png` land in `fret-pipeline/diagrams/`. Copy the figure you want to cite into `docs/figures/` so the thesis can pick it up. -2. **Plant change.** Edit `plant-model/pke_params.m` or dynamics in - `pke_th_rhs.m` → re-run `main.m` to verify behavior. Steady-state values - feed the `T_avg`, `P_crit` thresholds in the FRET requirements — if you - change them here, check the requirement predicates match. -3. **Verification work.** Continuous-mode verification (reachability, barrier - certificates) reads the guards from `fret-pipeline/specs/synthesis_config_v3.json` - and the dynamics from `plant-model/pke_th_rhs.m`. Results go in `reachability/`. +2. **Plant change.** Edit `code/src/pke_params.jl` or dynamics in + `code/src/pke_th_rhs.jl` → re-run `code/scripts/main_mode_sweep.jl` + to verify behavior. Steady-state values feed the `T_avg`, `P_crit` + thresholds via `reachability/predicates.json` — if you change + parameters here, the predicate concretization may need to follow. +3. **Verification work.** Continuous-mode verification reads the + FRET-side predicates from `reachability/predicates.json` and the + dynamics from `code/src/pke_th_rhs.jl`. Reach scripts in + `code/scripts/reach_*.jl` and `barrier_*.jl`. Results land back in + `reachability/*.mat` (gitignored). 4. **Thesis edit.** `cd thesis && && git commit && git push origin main`, then from the umbrella root: `git add thesis && git commit` to bump the submodule pointer. ## Where to look for what -- Reactor physics, steady-state, ODE dynamics → `plant-model/` +- Reactor physics, steady-state, ODE dynamics, controllers, reach, barrier → `code/` +- Predicate concretizations (the seam between FRET and reach) → `reachability/predicates.json` +- Reach analysis writeup → `reachability/WALKTHROUGH.md` - Controller synthesis, requirements, state machines → `fret-pipeline/` - Research narrative, math, claims, citations → `thesis/` - Shared figures (for the thesis or talks) → `docs/figures/` - How the layers compose → `docs/architecture.md` +- Project history / invention log → `journal/` +- Predicate-explorer GUI → `app/` ## Submodule rituals @@ -92,10 +100,13 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM). - FRET variables: `control_ = q_` for modes; everything else is an environment input. See `fret-pipeline/README.md` for the full naming rule. - Plant model units: internal SI (W, kg, °C); printed/plotted in °F. -- Continuous state vector: `y = [n; C1..C6; T_f; T_c; T_cold]`. Guards in - FRET should reference the algebraic outputs (`T_hot`, `T_avg`) or direct - states, not derived quantities. +- Continuous state vector: `x = [n; C1..C6; T_f; T_c; T_cold]`. Guards + in FRET should reference the algebraic outputs (`T_hot`, `T_avg`) or + direct states, not derived quantities. - `Q_sg(t)` is the bounded disturbance — never treat `T_cold` as an input. +- The continuous toolchain is **Julia**. The earlier MATLAB + implementation was ported on 2026-04-20 and deleted; recover via + git history if archaeologically needed. ## What's deliberately missing (so you don't go looking) @@ -105,5 +116,13 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM). diagram → hand-translate to Stateflow. Roadmapped. - Hardware bring-up (Ovation DCS) — scheduled for the integration thrust. `hardware/` is a placeholder. -- Continuous-mode verification — `reachability/` is empty. Choice of tool - (CORA / Flow* / SpaceEx / JuliaReach) still pending. +- Continuous-mode verification — partial. Linear reach for operation + mode discharges all 6 `inv2_holds` halfspaces (approximate, not + sound). Lyapunov-ellipsoid barrier fails structurally + (anisotropy ceiling — see `reachability/WALKTHROUGH.md`). Nonlinear + reach (TMJets) functional to 10 s; needs reduced-order PKE for + hours-long horizons. +- Polytopic / SOS barriers — flagged in WALKTHROUGH as the next-step + alternative to the failing quadratic Lyapunov barrier. +- Reduced-order (prompt-jump) PKE — known remedy for nonlinear reach + stiffness, not yet implemented. diff --git a/README.md b/README.md index 1566c24..83dff5d 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,23 @@ # pwr-hybrid-3-demo -Preliminary example for the HAHACS thesis — a verified hybrid controller for -a small modular PWR startup. Composes three layers into one demonstrable -pipeline: +Preliminary example for the HAHACS thesis — a verified hybrid controller +for a small modular PWR startup. Composes three layers into one +demonstrable pipeline: -- **Discrete layer** (`fret-pipeline/`): FRET natural-language requirements - → LTL → synthesized AIGER controller → state-machine diagram. -- **Continuous layer** (`plant-model/`): 10-state point kinetic equation + - thermal-hydraulics PWR model with bounded steam-generator heat removal as - the disturbance input. -- **Research context** (`thesis/`): the HAHACS PhD proposal that motivates - and formalizes the methodology. +- **Discrete layer** (`fret-pipeline/`): FRET natural-language + requirements → LTL → synthesized AIGER controller → state-machine + diagram. +- **Continuous layer** (`code/`): 10-state point kinetic equation + + thermal-hydraulics PWR model with bounded steam-generator heat + removal as the disturbance input. Controllers, linearization, LQR, + reach-tube propagator, Lyapunov barrier — all Julia. +- **Verification artifacts** (`reachability/`): predicate + concretizations (single source of truth in `predicates.json`) and + the standalone reach analysis writeup (`WALKTHROUGH.md`). +- **Research context** (`thesis/`): the HAHACS PhD proposal. +- **Lab journal** (`journal/`): chronological invention log in LaTeX. +- **Predicate explorer app** (`app/`): Pluto.jl notebook bridging + FRET predicates and continuous-state halfspaces. ## Layout @@ -18,14 +25,15 @@ pipeline: pwr-hybrid-3-demo/ CLAUDE.md AI-facing context and architecture map docs/ - architecture.md How the discrete and continuous layers compose + architecture.md How the layers compose figures/ Shared figures for thesis + talks fret-pipeline/ FRET → ltlsynt → AIGER → state machine - plant-model/ PWR point kinetics + thermal-hydraulics - reachability/ Continuous-mode verification (linear-model tube + Lyapunov barrier attempt; see README) - julia-port/ Parallel plant-model port + ReachabilityAnalysis.jl scaffold + code/ Plant model, controllers, reach (all Julia) + reachability/ predicates.json + WALKTHROUGH.md + app/ Pluto.jl predicate explorer + journal/ LaTeX lab notebook hardware/ Ovation HIL artifacts (TBD) - claude_memory/ Session notes by AI agents (distilled up into CLAUDE.md over time) + claude_memory/ Short AI-context notes thesis/ [submodule] PhD proposal presentations/ 2026DICE/ [submodule] DICE 2026 abstract @@ -50,37 +58,48 @@ python3 scripts/trace_aiger.py circuits/PWR_HYBRID_3_DRC.aag diagrams dot -Tpng diagrams/PWR_HYBRID_3_DRC_states.dot -o diagrams/PWR_HYBRID_3_DRC_states.png ``` -Run the plant model (MATLAB in `plant-model/` — Octave compatibility not tested since the LQR pieces landed): +Run the plant model and reach analysis: -```matlab -main % original single-scenario demo (null vs operation) -main_mode_sweep % all five DRC modes back-to-back, writes to ../docs/figures/ -test_linearize % Jacobian sanity check, saves linearization for reach +```bash +cd code +julia --project=. -e 'using Pkg; Pkg.instantiate()' # first time only + +julia --project=. scripts/main_mode_sweep.jl # all 5 DRC modes +julia --project=. scripts/reach_operation.jl # operation-mode linear reach +julia --project=. scripts/barrier_lyapunov.jl # Lyapunov barrier +julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier +julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup (10s cap) ``` -Run the reach artifacts (`reachability/`): +Open the predicate explorer: -```matlab -reach_operation % linear reach tube for operation-mode LQR -barrier_lyapunov % Lyapunov-ellipsoid barrier cert attempt (sweeps weights) +```bash +cd app +julia --project=. -e 'using Pluto; Pluto.run()' +# Browser opens; navigate to predicate_explorer.jl ``` -**Soundness note:** the current reach tube is the LINEAR model's tube; -it is not yet a sound over-approximation of the nonlinear plant. See -`reachability/README.md` § Soundness status. +**Soundness note:** the current reach tubes are over-approximations +of the LINEAR model, not sound over-approximations of the nonlinear +plant. See `reachability/README.md` and `reachability/WALKTHROUGH.md`. ## Prerequisites -- Python 3.10+ -- [Spot](https://spot.lre.epita.fr/) for `ltlsynt` (`brew install spot`) -- [Graphviz](https://graphviz.org/) for `dot` (`brew install graphviz`) -- MATLAB or GNU Octave for the plant model -- LaTeX (via `latexmk`) for the thesis submodule +- Julia 1.10+ (via `juliaup`). +- Python 3.10+ (FRET pipeline only). +- [Spot](https://spot.lre.epita.fr/) for `ltlsynt` (`brew install spot`). +- [Graphviz](https://graphviz.org/) for `dot` (`brew install graphviz`). +- LaTeX (via `latexmk`) for the thesis + journal builds. ## Further reading - `CLAUDE.md` — orientation for AI agents working in this repo - `docs/architecture.md` — how the layers compose +- `code/CLAUDE.md` — code architecture, conventions, validity range +- `code/README.md` — usage and dependencies +- `reachability/README.md` — reach scope, soundness status +- `reachability/WALKTHROUGH.md` — standalone analysis writeup +- `journal/README.md` — journal format conventions +- `journal/journal.tex` — the journal itself, dated entries - `thesis/CLAUDE.md` — the thesis project structure - `fret-pipeline/README.md` — FRET naming conventions and pipeline details -- `plant-model/README.md` — scenario setup and model equations diff --git a/julia-port/.gitignore b/code/.gitignore similarity index 100% rename from julia-port/.gitignore rename to code/.gitignore diff --git a/code/CLAUDE.md b/code/CLAUDE.md new file mode 100644 index 0000000..d9b1585 --- /dev/null +++ b/code/CLAUDE.md @@ -0,0 +1,185 @@ +# CLAUDE.md + +Guidance for Claude Code (and any AI agent) working in this subdirectory. + +See also: the parent repo's `CLAUDE.md` (living-documentation rule — +update this file when new knowledge is created), `claude_memory/` for +session notes, and `../journal/` for the long-form invention log. + +## What this is + +The full plant-model + reachability toolchain for the PWR_HYBRID_3 +preliminary example. All Julia — the earlier MATLAB implementation +was ported end-to-end on 2026-04-20 and deleted. A 10-state point +kinetic equation (PKE) model of a PWR coupled with lumped-parameter +thermal-hydraulics, plus mode-specific continuous controllers, linear +reach-set machinery, Lyapunov-barrier attempt, and scaffolding for +nonlinear reach via TMJets. + +The model is being developed for hybrid-systems reachability analysis. +The end goal is nonlinear reach (sound, not linearized) over per-mode +reach-avoid obligations, which will require a reduced-order PKE to +eliminate prompt-neutron stiffness. + +## Naming convention (thesis-aligned) + +Continuous state and control-theory notation matches the thesis: + +| Symbol | Meaning | +|---|---| +| `t` | time [s] | +| `x` | continuous state vector, `[n; C1..C6; T_f; T_c; T_cold]` | +| `plant` | parameter NamedTuple from `pke_params()` | +| `u` | control input: external rod reactivity [dk/k] | +| `Q_sg` | bounded disturbance: SG heat removal [W], as a callable `Q_sg(t)` | +| `ref` | reference/setpoint NamedTuple (e.g. `ref.T_avg`) | + +`T_avg = T_c` (algebraic alias); `T_hot = 2*T_c - T_cold` (linear +coolant profile assumption). + +## Running + +From the project root: + +```bash +cd code +julia --project=. -e 'using Pkg; Pkg.instantiate()' # first time only + +julia --project=. scripts/main_mode_sweep.jl # runs all 5 DRC modes +julia --project=. scripts/reach_operation.jl # operation-mode linear reach +julia --project=. scripts/barrier_lyapunov.jl # Lyapunov-ellipsoid barrier +julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier bounds +julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup reach (10s max, stiffness-limited) +julia --project=. scripts/sim_sanity.jl # operation-mode nominal sim +julia --project=. scripts/sim_heatup.jl # heatup nominal sim +``` + +Figures save to `../docs/figures/`. Reach results save to +`../reachability/*.mat` (gitignored). + +## Architecture + +``` +code/ + Project.toml Julia package manifest + README.md usage overview + src/ + pke_params.jl plant parameters and derived steady state + 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 + controllers/ + controllers.jl ctrl_null, shutdown, heatup, operation, + operation_lqr factory, scram + ctrl_heatup_unsat.jl saturation-disabled variant for reach + scripts/ + main_mode_sweep.jl runs every mode, saves all figures + sim_sanity.jl operation-mode Rodas5 sim (MATLAB-parity check) + sim_heatup.jl heatup Rodas5 sim + reach_operation.jl linear reach tube + inv2_holds check + barrier_lyapunov.jl Lyapunov barrier sweep + barrier_compare_OL_CL.jl OL vs CL barrier comparison + reach_heatup_nonlinear.jl TMJets nonlinear reach (10s horizon cap) +``` + +## Controller contract + +```julia +u = ctrl_mode(t, x, plant, ref) # all five args required +``` + +Pure function. Returns scalar `u` (external rod reactivity in `dk/k`). +All four signature args are required even if unused — lets the solver +swap controllers by function reference without caring which one it is. + +## Solver contract + +```julia +t, X, U = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan; x0=nothing, verbose=true) +``` + +Returns `(t::Vector, X::Matrix, U::Vector)` where `X[k, :]` is the +state at time `t[k]` and columns of X match the state vector ordering. +The underlying integrator is `Rodas5` from OrdinaryDiffEq (stiff-aware, +matches the MATLAB `ode15s` behavior we validated against). Optional +`x0` argument for starting from a non-equilibrium IC (heatup, shutdown). + +## Key design decisions + +- **`Q_sg` is the disturbance, never an actuator.** Cold leg + temperature is a dynamic state resulting from SG heat removal — the + reachability formulation needs `Q_sg ∈ [Q_min, Q_max]` as a bounded + disturbance, not a control input. +- **`u` is explicit.** `pke_th_rhs!` takes `u` directly instead of + reading it from the plant struct. +- **Controller gains live in the controller file.** Rough-out phase — + once we tune more than one mode we'll lift to a config struct. +- **`Rodas5` is required** because the system is stiff: + `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. +- **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 + struct-field accesses. Keep these in sync. + +## Conventions + +- One function per file where natural; group related functions in + one file when they're small (the controllers all live in one file). +- Internal model uses SI (W, kg, °C); printed/plotted temperatures + in °F. +- Reference setpoints in SI (so `ref.T_avg` is in °C). +- Controllers are pure functions — no persistent state. If a mode + needs integral action, augment `x`, don't hide state in globals. + (Exception: the LQR factory caches K in a closure, but that's + immutable after construction.) + +## Model validity range + +The feedback coefficients `alpha_f`, `alpha_c` are *linear* slopes +taken about the full-power operating point `(T_f0, T_c0)`. They become +unphysical when extrapolated far from that reference. True cold +shutdown (~50 °C everywhere) violates the linear-coefficient +assumption. + +**Trust range:** roughly ±50 °C around the operating point. +"Shutdown" in this demo means hot standby (~275 °C per +`predicates.json` `T_standby`), not cold shutdown. Extending to +cold-shutdown semantics requires piecewise-linear or table-lookup +temperature coefficients — separate modeling project. + +## For reachability work + +- `pke_th_rhs!(dx, x, t, plant, Q_sg, u)` is the `ẋ = f(x, u)` with + `Q_sg` as the disturbance `w`. +- Use `pke_linearize(plant; x_star, u_star, Q_star)` for `(A, B, B_w)` + at any operating point. +- Safe regions come from `../reachability/predicates.json` — three + layers: `operational_deadbands` (DRC transitions), + `safety_limits` (hard halfspaces), `mode_invariants` (conjunctions + used by per-mode reach). +- Do NOT try to verify the whole hybrid system at once. Pick one + mode, compute its reach set, discharge the per-mode obligation. +- Linear reach in `scripts/reach_operation.jl`; nonlinear (TMJets) + in `scripts/reach_heatup_nonlinear.jl`. + +## Robustness caveats (idealized in current artifacts) + +- **α_f, α_c are treated as known exactly.** The feedback-linearization + in `ctrl_heatup` assumes the controller's α matches the plant's; + real reactors have α drift (burnup, boron, xenon). Robust treatment + = α as bounded parametric uncertainty. +- **Saturation is a hybrid sub-mode.** `ctrl_heatup`'s `clamp(u, ...)` + is formally piecewise affine. Current reach treats it as dormant + (true for operation/LQR, near-true for the demo heatup trajectory). + A rigorous heatup reach has to model saturation regions explicitly. +- **Linear-model reach is not sound for the nonlinear plant.** The + linear-reach artifacts produce a sound over-approximation of the + LINEAR model's reach set, not of the plant's. To upgrade: + nonlinear reach directly (TMJets, once the reduced-order PKE + eliminates prompt-neutron stiffness), or Taylor-remainder inflation + on linear tubes. diff --git a/julia-port/Project.toml b/code/Project.toml similarity index 100% rename from julia-port/Project.toml rename to code/Project.toml diff --git a/code/README.md b/code/README.md new file mode 100644 index 0000000..b0ac634 --- /dev/null +++ b/code/README.md @@ -0,0 +1,59 @@ +# code + +Plant model, controllers, and reach-analysis toolchain for the +HAHACS preliminary example. All Julia. + +## What this is + +A 10-state coupled neutronics + thermal-hydraulics model (point kinetic +equations + lumped thermal loop) with continuous-mode controllers for +each of the DRC modes (shutdown, heatup, operation, scram), plus a +hand-rolled linear reach-tube propagator, a Lyapunov-ellipsoid barrier +attempt, and scaffolding for TMJets-based nonlinear reach. + +Ported from MATLAB on 2026-04-20 once the reach experiments made it +clear that Julia's stack (`OrdinaryDiffEq`, `MatrixEquations`, +`ReachabilityAnalysis`, `LazySets`, `@taylorize`) was the right tool +for everything we need going forward. The MATLAB originals are in +the git history. + +## Running + +First time: + +```bash +cd code +julia --project=. -e 'using Pkg; Pkg.instantiate()' +``` + +Subsequent: + +```bash +julia --project=. scripts/main_mode_sweep.jl # all 5 DRC modes, figures +julia --project=. scripts/reach_operation.jl # operation-mode linear reach +julia --project=. scripts/barrier_lyapunov.jl # Lyapunov barrier attempt +julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier +julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup (10s cap) +``` + +Figures save to `../docs/figures/`. Reach results save to +`../reachability/*.mat` (gitignored). + +## Structure + +See `CLAUDE.md` for the architectural overview and `../journal/` for +the invention-log-style narrative of how this code got written. + +## Dependencies + +From `Project.toml`: + +- `OrdinaryDiffEq` — ODE solver, Rodas5 for stiff systems. +- `MatrixEquations` — `arec` for LQR Riccati, `lyapc` for Lyapunov. +- `ReachabilityAnalysis` + `LazySets` — reach sets and set operations. +- `Plots` — figures (GR backend by default). +- `JSON` — read `../reachability/predicates.json`. +- `MAT` — save results. + +`Manifest.toml` is gitignored; regenerate locally on first +`Pkg.instantiate()`. diff --git a/julia-port/controllers/controllers.jl b/code/controllers/controllers.jl similarity index 100% rename from julia-port/controllers/controllers.jl rename to code/controllers/controllers.jl diff --git a/julia-port/controllers/ctrl_heatup_unsat.jl b/code/controllers/ctrl_heatup_unsat.jl similarity index 100% rename from julia-port/controllers/ctrl_heatup_unsat.jl rename to code/controllers/ctrl_heatup_unsat.jl diff --git a/julia-port/scripts/barrier_compare_OL_CL.jl b/code/scripts/barrier_compare_OL_CL.jl similarity index 100% rename from julia-port/scripts/barrier_compare_OL_CL.jl rename to code/scripts/barrier_compare_OL_CL.jl diff --git a/julia-port/scripts/barrier_lyapunov.jl b/code/scripts/barrier_lyapunov.jl similarity index 100% rename from julia-port/scripts/barrier_lyapunov.jl rename to code/scripts/barrier_lyapunov.jl diff --git a/julia-port/scripts/main_mode_sweep.jl b/code/scripts/main_mode_sweep.jl similarity index 100% rename from julia-port/scripts/main_mode_sweep.jl rename to code/scripts/main_mode_sweep.jl diff --git a/julia-port/scripts/reach_heatup_nonlinear.jl b/code/scripts/reach_heatup_nonlinear.jl similarity index 100% rename from julia-port/scripts/reach_heatup_nonlinear.jl rename to code/scripts/reach_heatup_nonlinear.jl diff --git a/julia-port/scripts/reach_operation.jl b/code/scripts/reach_operation.jl similarity index 100% rename from julia-port/scripts/reach_operation.jl rename to code/scripts/reach_operation.jl diff --git a/julia-port/scripts/sim_heatup.jl b/code/scripts/sim_heatup.jl similarity index 100% rename from julia-port/scripts/sim_heatup.jl rename to code/scripts/sim_heatup.jl diff --git a/julia-port/scripts/sim_sanity.jl b/code/scripts/sim_sanity.jl similarity index 100% rename from julia-port/scripts/sim_sanity.jl rename to code/scripts/sim_sanity.jl diff --git a/julia-port/src/load_predicates.jl b/code/src/load_predicates.jl similarity index 100% rename from julia-port/src/load_predicates.jl rename to code/src/load_predicates.jl diff --git a/julia-port/src/pke_linearize.jl b/code/src/pke_linearize.jl similarity index 100% rename from julia-port/src/pke_linearize.jl rename to code/src/pke_linearize.jl diff --git a/julia-port/src/pke_params.jl b/code/src/pke_params.jl similarity index 100% rename from julia-port/src/pke_params.jl rename to code/src/pke_params.jl diff --git a/julia-port/src/pke_solver.jl b/code/src/pke_solver.jl similarity index 100% rename from julia-port/src/pke_solver.jl rename to code/src/pke_solver.jl diff --git a/julia-port/src/pke_th_rhs.jl b/code/src/pke_th_rhs.jl similarity index 100% rename from julia-port/src/pke_th_rhs.jl rename to code/src/pke_th_rhs.jl diff --git a/julia-port/src/plot_pke_results.jl b/code/src/plot_pke_results.jl similarity index 100% rename from julia-port/src/plot_pke_results.jl rename to code/src/plot_pke_results.jl diff --git a/julia-port/src/reach_linear.jl b/code/src/reach_linear.jl similarity index 100% rename from julia-port/src/reach_linear.jl rename to code/src/reach_linear.jl diff --git a/journal/journal.tex b/journal/journal.tex index 4704e4f..6a19ba0 100644 --- a/journal/journal.tex +++ b/journal/journal.tex @@ -33,6 +33,15 @@ Each section is a dated session. Sessions are written in two styles: content that should be expanded in a later A-pass. \end{itemize} +\textbf{File-path archaeology note.} Entries dated on or before +2026-04-20 (afternoon) refer to MATLAB files under \texttt{plant-model/} +and a MATLAB-Julia split where Julia code lived under \texttt{julia-port/}. +The 2026-04-20 evening mega-session ported everything to Julia and +deleted MATLAB; the result is the unified \texttt{code/} directory. +Path references in earlier entries are historically accurate for their +date; recover the corresponding sources via \texttt{git log} or the +parent commit at the time. + Callout boxes signal specific content types: \begin{derivation} diff --git a/julia-port/README.md b/julia-port/README.md deleted file mode 100644 index 55e46a4..0000000 --- a/julia-port/README.md +++ /dev/null @@ -1,40 +0,0 @@ -# Julia port - -Parallel implementation of the plant model (`../plant-model/`) in Julia, -intended as an eventual landing spot for reachability if we outgrow the -MATLAB / hand-rolled tooling. - -## Status - -- `src/pke_params.jl`, `src/pke_th_rhs.jl`, `src/pke_linearize.jl` — - functional, match MATLAB term-for-term. -- `controllers/controllers.jl` — all five modes ported - (null, shutdown, heatup, operation, scram). LQR factory included via - MatrixEquations.jl. -- `scripts/sim_sanity.jl` — closed-loop simulation matches MATLAB to 3 - decimals on the 100% → 80% `Q_sg` step. ✅ passing. -- `scripts/reach_operation.jl` — ❌ stub. ReachabilityAnalysis.jl - algorithms blow up on this stiff, badly-conditioned system. See the - file header for diagnosis and planned fix (state rescaling). - -## When to prefer Julia over MATLAB - -Today: nowhere. The sanity path exists so we don't regret the eventual -port. - -Once the reach scaling is resolved: -- Nonlinear reach for the P controller in operation mode (CORA / - JuliaReach territory where MATLAB's linearization doesn't suffice). -- Heatup reach with its time-varying reference. -- Parametric studies where MATLAB license fees / CI friction matter. - -## Running - -```bash -cd julia-port -julia --project=. -e 'using Pkg; Pkg.instantiate()' -julia --project=. scripts/sim_sanity.jl -``` - -First run will precompile the dependency stack (OrdinaryDiffEq, -ReachabilityAnalysis, LazySets — a few minutes). diff --git a/plant-model/CLAUDE.md b/plant-model/CLAUDE.md deleted file mode 100644 index a00945f..0000000 --- a/plant-model/CLAUDE.md +++ /dev/null @@ -1,196 +0,0 @@ -# CLAUDE.md - -Guidance for Claude Code (and any AI agent) working in this subdirectory. - -See also: the parent repo's `CLAUDE.md` (living-documentation rule — update -this file when new knowledge is created) and `claude_memory/` for session -notes that haven't yet earned a place here. - -## What this is - -A 10-state point kinetic equation (PKE) model of a PWR coupled with -lumped-parameter thermal-hydraulics, plus mode-specific continuous -controllers. Implements the *plant* and the *tactical* layer of the HAHACS -methodology — see `../thesis/3-research-approach/approach.tex` sections on -transitory, stabilizing, and expulsory modes. - -The model is being developed for hybrid-systems reachability analysis. The -end goal is barrier-certificate and reachability methods (in `../reachability/`) -to prove each continuous mode satisfies its obligations under bounded -`Q_sg(t)` variations. - -## Naming convention (thesis-aligned) - -Continuous state and control-theory notation matches the thesis: - -| Symbol | Meaning | -|---|---| -| `t` | time [s] | -| `x` | continuous state vector, `[n; C1..C6; T_f; T_c; T_cold]` | -| `plant` | parameter struct from `pke_params()` (was `p` — renamed to avoid collision with `p_i` predicates) | -| `u` | control input: external rod reactivity [dk/k] | -| `Q_sg` | bounded disturbance: SG heat removal [W], as a function handle `Q_sg(t)` | -| `ref` | reference/setpoint struct (e.g. `ref.T_avg`) | - -`T_avg = T_c` (algebraic alias); `T_hot = 2*T_c - T_cold` (linear coolant -profile assumption). - -## Running - -From inside MATLAB: - -```matlab -main % runs the default closed-loop scenario -``` - -From a terminal (figures appear as their own windows): - -```bash -matlab -nodesktop -nosplash -r "main" # drops at >> prompt; `exit` to quit -matlab -r "main" & # full IDE, backgrounded -matlab -batch "main" # headless sanity check (no figures shown) -``` - -`main.m` configures the scenario (plant, `Q_sg`, ref, tspan), picks a -controller, and calls `pke_solver`. Tested on MATLAB R2025b; works on any -MATLAB with `ode15s` (R2020b or newer). - -## Architecture - -``` -plant-model/ - pke_params.m plant parameters and derived steady state - pke_initial_conditions.m analytic equilibrium x0 from plant - pke_th_rhs.m dynamics: dx/dt = f(t, x, plant, Q_sg, u) - pke_linearize.m numerical (A, B, B_w) at any x_star - pke_solver.m closed-loop driver (optional x0 arg) - plot_pke_results.m 4-panel results plot - load_profile.m canned Q_sg demand shapes - main.m entry point — single-mode scenarios - main_mode_sweep.m runs all DRC modes back-to-back - test_linearize.m sanity-check Jacobians vs nonlinear sim - controllers/ - ctrl_null.m u = 0 (baseline, no feedback) - ctrl_shutdown.m hot standby: u = -5*beta (constant) - ctrl_heatup.m ramp T_avg: feedback lin. + P on ramp, saturated u - ctrl_operation.m stabilizing: P on T_avg - ctrl_operation_lqr.m full-state LQR around x_op (persistent-cached K) - ctrl_scram.m emergency: u = -8*beta (constant) -``` - -**Controller contract:** - -```matlab -u = ctrl_(t, x, plant, ref) -``` - -Pure function. Returns scalar `u` (external rod reactivity in `dk/k`). All -four signature args are required even if unused — lets the solver swap -controllers by function handle without caring which one it is. - -**Solver contract:** - -```matlab -[t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan) -``` - -Builds the closed-loop RHS internally, calls `ode15s`, then reconstructs the -control trajectory `U` by re-evaluating `ctrl_fn` at each returned time step. -For event-driven or hysteretic controllers, this reconstruction may miss -interior solver evaluations — revisit if / when that matters. - -## Key design decisions - -- **`Q_sg` is the disturbance, never an actuator.** Cold leg temperature is - a dynamic state resulting from SG heat removal — the reachability - formulation needs `Q_sg ∈ [Q_min, Q_max]` as a bounded disturbance, not a - control input. -- **`u` is explicit.** `pke_th_rhs` takes `u` directly instead of reading - `p.rho_ext(t)` from the params struct. This decouples plant from - controller and matches the thesis's `f(x, u)` formulation. -- **Controller gains live in the controller file.** Rough-out phase — we'll - lift to a config struct once we're tuning more than one mode. The gain - in `ctrl_operation` (`Kp = 1e-4 /K`) is chosen to roughly double the - effective moderator coefficient, not for precision tracking. -- **`ode15s` is required** because the system is stiff: `Lambda` (~10⁻⁴ s) - vs thermal time constants (~10–100 s). - -## Conventions - -- One function per file, each file has a single responsibility. -- Internal model uses SI (W, kg, °C); all printed/plotted temperatures in °F. -- Reference setpoints in SI (so `ref.T_avg` is in °C). -- Controllers are pure functions — no persistent state. If a mode ever - needs integral action, augment `x` rather than hiding state in globals. - -## Model validity range - -The feedback coefficients `alpha_f`, `alpha_c` are *linear* slopes taken -about the full-power operating point `(T_f0, T_c0)`. They become -unphysical when extrapolated far from that reference — e.g. at true -cold shutdown (~50 C everywhere), the model predicts ~+5 $ intrinsic -reactivity that real reactors do not exhibit because temperature -coefficients saturate and flip sign in the cold-unborated state. - -**Trust range:** roughly ±50 C around the operating point. "Shutdown" -in this demo means hot standby (~290 C, n≈0), not cold shutdown. -Extending to cold-shutdown semantics requires piecewise-linear or -table-lookup temperature coefficients and is out of scope for the -running example. - -## For reachability work - -- `pke_th_rhs(t, x, plant, Q_sg, u)` is the `ẋ = f(x, u)` with `Q_sg` as the - disturbance `w`. For a given mode, substitute `u = ctrl_(t, x, plant, ref)` - to get the closed-loop `ẋ = f_cl(x, w)`. -- Use `pke_linearize(plant, x_star, u_star, Q_star)` for `(A, B, B_w)` at - any operating point; returns central-finite-difference Jacobians. The - result at the full-power point is saved to - `../reachability/linearization_at_op.mat` by `test_linearize.m`. -- Safe regions come from the FRET spec's guards (see - `../fret-pipeline/specs/synthesis_config_v3.json`) — these define - `X_entry`, `X_exit`, `X_safe` per the thesis approach. Numerical - predicate thresholds still need to be pinned (see - `../claude_memory/` notes for the discussion). -- Do NOT try to verify the whole hybrid system at once. Pick one mode, - compute its reach set, discharge the per-mode obligation. The - compositionality argument in the thesis is what makes this tractable. -- First reach set is in `../reachability/reach_operation.m` (operation - mode under LQR); a Lyapunov-ellipsoid barrier-cert attempt is in - `../reachability/barrier_lyapunov.m`. - -## For control design - -- The LQR gain in `ctrl_operation_lqr.m` is cached in a persistent - variable. If you change `Q_lqr` or `R_lqr` or mutate `plant`, restart - MATLAB to clear the cache. -- `ctrl_heatup.m` uses feedback linearization + saturated P. No - integrator on purpose — adding PI would double the plant's intrinsic - thermal-mass integration and make anti-windup a hybrid mode, both - bad for reach. -- Starting heatup from a truly zero-power IC produces a large lag + - power-spike overshoot in the simulation. Physical resolution: the - reactor should be taken critical in shutdown mode (at ~0.1% power) - before DRC transitions to heatup. `main_mode_sweep.m` uses this IC. - -## Robustness caveats (idealized in current artifacts) - -- **α_f, α_c are treated as known exactly.** In reality α_f drifts - ~20% over burnup; α_c spans ~10x across soluble-boron dilution over - a cycle; xenon adds 2-3 $ reactivity on its own timescale. The - feedback-linearization in `ctrl_heatup.m` assumes the controller's - α matches the plant's; if not, the clean `rho_total = Kp*e` - property degrades to `Kp*e + delta*alpha*dT`, and the P term must - absorb the residual. Stabilization still holds but reach analysis - should eventually treat α as a bounded parametric uncertainty. -- **Saturation is a hybrid sub-mode.** The `sat(u, u_min, u_max)` in - `ctrl_heatup.m` is formally piecewise affine. Current reach - treats it as dormant (true for operation/LQR, near-true for the - demo heatup trajectory). A rigorous heatup reach has to model - the saturation regions explicitly. -- **Linear-model reach is not sound for the nonlinear plant.** The - reach artifacts in `../reachability/` use the linearization; the - result is a sound over-approximation of the LINEAR model's reach, - not of the plant's. To upgrade: nonlinear reach directly, or - linear reach + Taylor-remainder inflation. See - `../reachability/README.md` § Soundness status. diff --git a/plant-model/README.md b/plant-model/README.md deleted file mode 100644 index 7d2901a..0000000 --- a/plant-model/README.md +++ /dev/null @@ -1,92 +0,0 @@ -# plant-model - -PWR plant model (point kinetics + lumped thermal-hydraulics) and -mode-specific continuous controllers for the HAHACS preliminary example. - -## Overview - -A 10-state coupled neutronics + thermal-hydraulics model in MATLAB: - -- 6 delayed neutron precursor groups (U-235 thermal fission, Keepin) -- Lumped fuel, core coolant, and SG/cold-leg thermal nodes -- Steam generator heat removal `Q_sg(t)` as the bounded disturbance input -- Doppler and moderator temperature reactivity feedback -- External rod reactivity `u` as the controllable input - -State vector: `x = [n; C1..C6; T_f; T_c; T_cold]` (10 states). See -`CLAUDE.md` for the naming convention. - -## Quick Start - -Open MATLAB in this directory and run: - -```matlab -main -``` - -The default scenario runs two simulations of a 100% → 80% SG demand step: -once with `ctrl_null` (plant feedback only) and once with `ctrl_operation` -(proportional rod reactivity on T_avg error), and plots the comparison. - -## Files - -| File | Role | -|------|------| -| `main.m` | Entry point — scenario config and run | -| `pke_params.m` | Plant parameters and steady-state derivation | -| `pke_th_rhs.m` | Dynamics `ẋ = f(t, x, plant, Q_sg, u)` | -| `pke_initial_conditions.m` | Analytic steady-state `x0` | -| `pke_solver.m` | Closed-loop driver — takes a controller function handle | -| `plot_pke_results.m` | 4-panel results plot | -| `load_profile.m` | SG heat demand shapes | -| `controllers/ctrl_null.m` | `u = 0` baseline | -| `controllers/ctrl_operation.m` | Stabilizing mode: P on `T_avg` | - -## Controllers - -Controllers share a single signature: - -```matlab -u = ctrl_(t, x, plant, ref) -``` - -Returns scalar `u` (external rod reactivity in `dk/k`). The solver swaps -controllers via function handle: - -```matlab -[t, X, U] = pke_solver(plant, Q_sg, @ctrl_operation, ref, tspan); -``` - -Additional modes (`ctrl_heatup`, `ctrl_scram`, `ctrl_shutdown`) will land in -`controllers/` following the same signature. - -## Running Different Scenarios - -Swap `Q_sg` in `main.m`: - -```matlab -% Step down to 90% at t = 10s -Q_sg = @(t) plant.P0 * (1.0 - 0.1 * (t >= 10)); - -% Interpolated time series -t_data = [0, 100, 200, 300]; -q_data = [1.0, 0.85, 0.9, 1.0] * plant.P0; -Q_sg = @(t) interp1(t_data, q_data, t, 'linear', 'extrap'); -``` - -Swap the controller: - -```matlab -[t, X, U] = pke_solver(plant, Q_sg, @ctrl_null, [], tspan); -``` - -Change the reference (for modes that use one): - -```matlab -ref.T_avg = plant.T_c0 + 5; % track 5 C above nominal -``` - -## Requirements - -MATLAB (R2020b or newer, tested on R2025b). Uses `ode15s` from base MATLAB -— no toolboxes required. diff --git a/plant-model/context.md b/plant-model/context.md deleted file mode 100644 index 8fcef54..0000000 --- a/plant-model/context.md +++ /dev/null @@ -1,88 +0,0 @@ -# PKE Playground — LLM Context - -This document provides context for an LLM working on this codebase. It explains the model, design decisions, and where the project is headed. - -## What This Is - -A point kinetic equation (PKE) model coupled with lumped-parameter thermal-hydraulics, implemented in MATLAB/Octave. The model represents a simplified PWR (pressurized water reactor) core and primary loop. - -## Why It Exists - -This is a workable plant model being developed for **hybrid systems reachability analysis**. The end goal is to use barrier certificate methods to prove that under bounded variations in steam generator (SG) heat removal, the reactor stays within a defined safe operating region (bounded temperatures, power levels, etc.). - -The model needs to be: -- Physically motivated (not just a toy system) -- Clean and state-space oriented (explicit state vector, clear inputs/outputs) -- Structured so that `Q_sg(t)` is a bounded disturbance input for reachability tools - -## The Model - -### State Vector (10 states) - -``` -y = [n; C1; C2; C3; C4; C5; C6; T_f; T_c; T_cold] -``` - -- `n` — normalized neutron population (1.0 = nominal) -- `C1..C6` — six delayed neutron precursor group concentrations (Keepin U-235 thermal fission data) -- `T_f` — lumped fuel temperature [C] -- `T_c` — average core coolant temperature [C] -- `T_cold` — cold leg / SG outlet coolant temperature [C] - -### Algebraic Outputs - -- `T_hot = 2*T_c - T_cold` (linear coolant temperature profile through the core) -- `T_avg = T_c` - -### Equations - -**Neutronics:** -``` -dn/dt = (rho - beta) / Lambda * n + sum(lambda_i * C_i) -dC_i/dt = beta_i / Lambda * n - lambda_i * C_i -``` - -**Thermal-hydraulics (three lumped nodes):** -``` -M_f*c_f * dT_f/dt = P0*n - hA*(T_f - T_c) [fuel] -M_c*c_c * dT_c/dt = hA*(T_f - T_c) - 2*W*c_c*(T_c - T_cold) [core coolant] -M_sg*c_c * dT_cold/dt = W*c_c*(T_hot - T_cold) - Q_sg(t) [SG/cold leg] -``` - -**Reactivity feedback:** -``` -rho = rho_ext(t) + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0) -``` - -### Why Q_sg Is the Input (Not T_cold) - -In a real PWR, the steam generator removes heat from the primary loop based on secondary-side demand. The cold leg temperature is a *result* of that heat removal, not a boundary condition you get to set. Treating `Q_sg(t)` as the external input and letting `T_cold` be a dynamic state is physically correct and maps naturally to the reachability problem: "given that SG demand is somewhere in [Q_min, Q_max], what states can the reactor reach?" - -### Why T_hot = 2*T_c - T_cold - -This comes from the linear coolant temperature profile assumption through the core. If coolant enters at `T_cold` and exits at `T_hot`, and the average is `T_c`, then `T_c = (T_hot + T_cold) / 2`, so `T_hot = 2*T_c - T_cold`. This avoids needing a separate hot leg state while still tracking both temperatures. - -### Why ode15s - -The PKE system is stiff. The prompt neutron generation time `Lambda` (~10^-4 s) is orders of magnitude smaller than the thermal time constants (~10-100 s) and precursor decay constants (~0.01-3 s^-1). An explicit solver would need extremely small timesteps; `ode15s` handles this naturally. - -## File Structure - -Each file has one purpose: - -| File | Purpose | -|------|---------| -| `pke_solver.m` | Main driver script — scenario config, solve, print | -| `pke_params.m` | All plant parameters and steady-state derivation | -| `pke_th_rhs.m` | ODE right-hand side `f(t, y)` — this is the dynamics | -| `pke_initial_conditions.m` | Analytic steady-state initial condition | -| `load_profile.m` | SG heat demand time series (swappable) | -| `plot_pke_results.m` | 4-panel results plotting | - -## For Future Reachability Work - -- `pke_th_rhs.m` is your `dx/dt = f(x, w)` where `w = Q_sg` -- Linearize around the steady state from `pke_params.m` to get `A`, `B` matrices -- The safe region will be a polytope or sublevel set on `[n, T_f, T_hot, T_cold]` (or a subset) -- `Q_sg` is the bounded disturbance: `Q_sg in [Q_min, Q_max]` around `P0` -- Barrier certificate `B(x)` would certify that trajectories starting in the safe set stay there for all admissible `Q_sg(t)` diff --git a/plant-model/controllers/ctrl_heatup.m b/plant-model/controllers/ctrl_heatup.m deleted file mode 100644 index 74371a3..0000000 --- a/plant-model/controllers/ctrl_heatup.m +++ /dev/null @@ -1,82 +0,0 @@ -function u = ctrl_heatup(t, x, plant, ref) -% CTRL_HEATUP Ramp T_avg toward a target at a bounded rate. -% -% Structure: -% u_ff = -alpha_f*(T_f - T_f0) - alpha_c*(T_c - T_c0) % cancel feedback -% T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target) % ramped reference -% u_unsat = u_ff + Kp * (T_ref - T_avg) % P on error -% u = sat(u_unsat, ref.u_min, ref.u_max) % bounded rod worth -% -% Why saturation: -% Without it the P gain can push u toward prompt-supercritical as the -% cold-hot feedback bias unwinds late in the ramp. Capping u at -% +0.5*beta guarantees rho_total < beta (below prompt), which in -% turn bounds the neutron-kinetics excursion rate for reachability. -% -% IMPORTANT for reach analysis: -% The sat() is a 3-mode piecewise-affine system (sat-low / linear / -% sat-high). Under linear reach assumptions it must either be -% (a) proven dormant (u_unsat stays in [u_min, u_max] across the -% reach set — trivial to check, expensive to over-approximate -% tightly), or -% (b) handled as an explicit hybrid automaton nested under the DRC -% mode, with transitions when u_unsat crosses the saturation -% bounds. -% The current reach_operation.m assumes (a) implicitly. Heatup -% reach will need option (b) because u_unsat is near the +0.5*beta -% bound during early ramp. -% -% Why no integrator: -% Ramp tracking has a structural lag proportional to ramp_rate / Kp_eff. -% Acceptable because the DRC exits heatup on a predicate window -% (t_avg_in_range & p_above_crit), not on zero steady-state error. -% Adding PI would double-count the intrinsic plant integrator -% (thermal mass) and make anti-windup a hybrid transition. -% -% IMPORTANT caveat on the cancellation u_ff: -% The feedback linearization works only if the controller's values -% of alpha_f, alpha_c match the plant's. In reality alpha drifts: -% alpha_f ~20% over burnup, alpha_c by ~10x across boron dilution, -% plus xenon. With alpha_true = alpha_nom*(1+delta), the -% cancellation leaves residual reactivity delta*alpha*dT that the -% P term cleans up — stabilization still holds, but the clean -% "rho_total = Kp*e" property is gone. A robust deployment would -% treat alpha as an interval and propagate parametric uncertainty -% through reach (zonotope with parameter generators), or add -% adaptive alpha estimation. Idealized here. -% -% Inputs: -% t - time [s] -% x - state vector (10 x 1) -% plant - parameter struct (alpha_f, alpha_c, T_f0, T_c0 used) -% ref - struct with fields: -% .T_start starting T_avg [C] -% .T_target final T_avg [C] -% .ramp_rate desired dT_avg/dt [C/s] -% .u_min (optional) lower saturation [dk/k]; default -5*beta -% .u_max (optional) upper saturation [dk/k]; default +0.5*beta - - Kp = 1e-4; % [dk/k per K] - - T_f = x(8); - T_c = x(9); - T_avg = T_c; - - % Feedforward: cancel intrinsic temperature feedback so rho_total = Kp*e - % (before saturation). - u_ff = -plant.alpha_f * (T_f - plant.T_f0) ... - -plant.alpha_c * (T_avg - plant.T_c0); - - % Ramped reference, clamped at target. - T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target); - - e = T_ref - T_avg; - u_unsat = u_ff + Kp * e; - - % Saturation bounds (defaults keep rod worth subcritical-prompt). - if isfield(ref, 'u_min'), u_min = ref.u_min; else, u_min = -5 * plant.beta; end - if isfield(ref, 'u_max'), u_max = ref.u_max; else, u_max = 0.5 * plant.beta; end - - u = min(max(u_unsat, u_min), u_max); - -end diff --git a/plant-model/controllers/ctrl_null.m b/plant-model/controllers/ctrl_null.m deleted file mode 100644 index 8f774e3..0000000 --- a/plant-model/controllers/ctrl_null.m +++ /dev/null @@ -1,13 +0,0 @@ -function u = ctrl_null(t, x, plant, ref) -% CTRL_NULL No-op controller. Returns u = 0. -% -% Use for unforced-plant simulations (feedback-only response to Q_sg) or -% as a baseline to compare against closed-loop modes. -% -% Signature matches all other controllers: -% u = ctrl_(t, x, plant, ref) -% so it can be dropped into pke_solver as a function handle. - - u = 0; - -end diff --git a/plant-model/controllers/ctrl_operation.m b/plant-model/controllers/ctrl_operation.m deleted file mode 100644 index 2090bb8..0000000 --- a/plant-model/controllers/ctrl_operation.m +++ /dev/null @@ -1,31 +0,0 @@ -function u = ctrl_operation(t, x, plant, ref) -% CTRL_OPERATION Stabilizing mode controller: hold T_avg at ref.T_avg. -% -% Proportional feedback on external reactivity (rod worth) driven by -% temperature error: -% u = Kp * (ref.T_avg - T_avg) -% -% Sign check: when T_avg > ref.T_avg, error < 0, u < 0 → rods insert, -% power drops, T_avg falls. Correct polarity for a PWR. -% -% Gain Kp is a rough-out value. Scale reasoning: -% - alpha_c ~ -1e-4 /K means the moderator already supplies 1e-4/K of -% restoring reactivity per Kelvin of temperature error. -% - Kp = 1e-4 /K roughly doubles the effective moderator coefficient — -% it damps transients without overwhelming the intrinsic feedback. -% - Real PWR rod controllers use much smaller gains with a deadband. -% We can retune or switch to PI once we care about steady-state offset. -% -% Inputs: -% t - time [s] (unused in this controller, kept for signature) -% x - state vector (10 x 1) -% plant - parameter struct (unused here; reserved for future gain-scheduling) -% ref - struct with field .T_avg = target average coolant temperature [C] - - Kp = 1e-4; % [dk/k per K] - - T_avg = x(9); - e = ref.T_avg - T_avg; - u = Kp * e; - -end diff --git a/plant-model/controllers/ctrl_operation_lqr.m b/plant-model/controllers/ctrl_operation_lqr.m deleted file mode 100644 index cd90691..0000000 --- a/plant-model/controllers/ctrl_operation_lqr.m +++ /dev/null @@ -1,53 +0,0 @@ -function u = ctrl_operation_lqr(t, x, plant, ref) -% CTRL_OPERATION_LQR Full-state feedback for operation mode. -% -% u = -K * (x - x_op), K from LQR on (A, B) at x_op. -% -% Advantages over ctrl_operation (plain P on T_avg): -% - Damps the fast neutronics and the slow thermal response jointly. -% - Adds NO state (K is 1x10, feedback is memoryless). -% - Closed-loop A_cl = A - B*K is trivially propagable for reachability. -% -% Q/R are persistently cached on the first call. If the caller mutates -% plant after the first invocation, gains become stale — clear `persistent` -% by restarting MATLAB. That's a reach-analysis bookkeeping issue, not a -% runtime one. -% -% Weighting rationale: -% Q diagonal: -% n -> 1 (power tracking matters) -% C1..C6 -> 1e-3 (precursor deviations are informational) -% T_f -> 1e-2 (care about it but not as much as coolant) -% T_c -> 1e2 (primary objective) -% T_cold -> 1 (secondary but consequential) -% R = 1e6 (discourage large rod movements — one Kp of |u| ~ 1e-3) -% -% Inputs: -% t, x, plant, ref - standard signature. ref is unused; the setpoint -% is baked in as x_op. If you want to track a -% different operating point, regenerate K there. - - persistent K x_op_cached - - if isempty(K) - x_op_cached = pke_initial_conditions(plant); - [A, B, ~] = pke_linearize(plant, x_op_cached, 0, plant.P0); - - Q = diag([1, ... - 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, ... - 1e-2, 1e2, 1]); - R = 1e6; - - % MATLAB's lqr requires Control System Toolbox. Fall back to - % icare (in base MATLAB from R2019a) if lqr is unavailable. - try - K = lqr(A, B, Q, R); - catch - [~, ~, K_ric] = icare(A, B, Q, R); - K = K_ric; - end - end - - u = -K * (x - x_op_cached); - -end diff --git a/plant-model/controllers/ctrl_scram.m b/plant-model/controllers/ctrl_scram.m deleted file mode 100644 index b07b72a..0000000 --- a/plant-model/controllers/ctrl_scram.m +++ /dev/null @@ -1,20 +0,0 @@ -function u = ctrl_scram(t, x, plant, ref) -% CTRL_SCRAM Emergency shutdown: max negative rod worth inserted. -% -% Applies the full scram rod worth instantaneously. In a real plant the -% rods free-fall in ~2-3 seconds; this idealized version steps to the -% final worth at t = t_scram. Good enough for first-pass reachability, -% conservative for safety arguments (faster insertion => less excursion). -% -% u = -8 * beta (~8 $, in the typical 8-10 $ regulatory band) -% -% Same feedback-linearization comment as ctrl_shutdown: this dominates -% any plausible feedback the reactor can produce en route from any -% power/temperature state. -% -% Inputs: -% t, x, plant, ref - standard signature; all unused here. - - u = -8 * plant.beta; - -end diff --git a/plant-model/controllers/ctrl_shutdown.m b/plant-model/controllers/ctrl_shutdown.m deleted file mode 100644 index a56014b..0000000 --- a/plant-model/controllers/ctrl_shutdown.m +++ /dev/null @@ -1,21 +0,0 @@ -function u = ctrl_shutdown(t, x, plant, ref) -% CTRL_SHUTDOWN Hot-standby / cold-shutdown passive controller. -% -% Holds the reactor deeply subcritical with rods fully in. No feedback — -% just a constant negative external reactivity that dominates any -% temperature-feedback contribution the reactor could plausibly produce -% in its accessible state set. -% -% u = -5 * beta (~5 $ subcritical) -% -% Sign check: rho_total = u + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0). -% At cold IC (T < T0) the feedback contribution is positive (roughly -% +2.8 $ at T = T_cold0 everywhere), so -5 $ leaves the reactor at about -% -2.2 $ — comfortably subcritical. -% -% Inputs: -% t, x, plant, ref - standard signature; all unused here. - - u = -5 * plant.beta; - -end diff --git a/plant-model/load_profile.m b/plant-model/load_profile.m deleted file mode 100644 index 989bebd..0000000 --- a/plant-model/load_profile.m +++ /dev/null @@ -1,26 +0,0 @@ -function f = load_profile(t) -% LOAD_PROFILE Returns fractional SG heat demand at time t. -% -% Output is in [0, 1] where 1.0 = 100% of nominal power P0. -% Multiply by P0 to get Q_sg in watts: Q_sg = @(t) P0 * load_profile(t) -% -% Default profile (load-following): -% 0 - 30s: hold at 100% -% 30 - 130s: ramp 100% -> 80% -% 130 - 350s: hold at 80% -% 350 - 450s: ramp 80% -> 100% -% 450+: hold at 100% - - if t < 30 - f = 1.0; - elseif t < 130 - f = 1.0 - 0.2 * (t - 30) / 100; - elseif t < 350 - f = 0.8; - elseif t < 450 - f = 0.8 + 0.2 * (t - 350) / 100; - else - f = 1.0; - end - -end diff --git a/plant-model/main.m b/plant-model/main.m deleted file mode 100644 index 113473a..0000000 --- a/plant-model/main.m +++ /dev/null @@ -1,50 +0,0 @@ -%% main.m — closed-loop entry point -% -% Scenario: a 100% → 80% SG demand step. Compare the plant's intrinsic -% feedback response (ctrl_null) against the stabilizing operation-mode -% controller (ctrl_operation) holding T_avg at the nominal setpoint. -% -% Controller signature (thesis naming): -% u = ctrl_(t, x, plant, ref) -% where x is the continuous state vector and u is the external rod reactivity. - -clear; clc; close all; -addpath('controllers'); - -%% ===== Plant ===== -plant = pke_params(); - -%% ===== Scenario ===== - -% SG heat removal: step down to 80% at t = 30s, hold. -Q_sg = @(t) plant.P0 * (1.0 - 0.2 * (t >= 30)); - -% Reference for the operation-mode controller. -ref = struct(); -ref.T_avg = plant.T_c0; % hold T_avg at steady-state value - -% Simulation window. -tspan = [0, 600]; - -%% ===== Baseline: unforced feedback response ===== -fprintf('\n----- Run 1: ctrl_null (plant feedback only) -----\n'); -[t1, X1, U1] = pke_solver(plant, Q_sg, @ctrl_null, [], tspan); - -%% ===== Closed loop: operation mode ===== -fprintf('\n----- Run 2: ctrl_operation (P on T_avg) -----\n'); -[t2, X2, U2] = pke_solver(plant, Q_sg, @ctrl_operation, ref, tspan); - -%% ===== Plots ===== -plot_pke_results(t1, X1, U1, plant, Q_sg, 'ctrl_null (baseline)'); -plot_pke_results(t2, X2, U2, plant, Q_sg, 'ctrl_operation (P on T_avg)'); - -%% ===== Quick comparison figure ===== -CtoF = @(T) T*9/5 + 32; -figure('Position', [100 50 1000 450], 'Name', 'T_avg comparison'); -plot(t1, CtoF(X1(:,9)), 'b-', 'LineWidth', 1.5); hold on; -plot(t2, CtoF(X2(:,9)), 'r-', 'LineWidth', 1.5); -yline(CtoF(plant.T_c0), 'k--', 'LineWidth', 1.0); -xlabel('Time [s]'); ylabel('T_{avg} [F]'); -legend('ctrl\_null', 'ctrl\_operation', 'setpoint', 'Location', 'best'); -title('T_{avg} tracking: baseline vs stabilizing controller'); -grid on; diff --git a/plant-model/main_mode_sweep.m b/plant-model/main_mode_sweep.m deleted file mode 100644 index bbe33c0..0000000 --- a/plant-model/main_mode_sweep.m +++ /dev/null @@ -1,138 +0,0 @@ -%% main_mode_sweep.m — demo each DRC continuous-mode controller -% -% One scenario per mode, from an initial condition plausible for that -% mode. Figures saved to ../docs/figures/mode_sweep_*.png so the runs -% are reviewable without the MATLAB IDE. -% -% Modes, in the order DRC visits them starting from shutdown: -% 1. shutdown — deep subcritical, cold IC, Q_sg = 0 -% 2. heatup — same cold IC, ramped T_avg reference -% 3. operation — operating IC, 100% -> 80% Q_sg step -% 4. scram — operating IC, rods slammed in, Q_sg decays to 0 -% -% Each run produces plot_pke_results 4-panel plus saves a lightweight -% summary figure. - -clear; clc; close all; -addpath('controllers'); -addpath(fullfile('..', 'reachability')); % for load_predicates - -plant = pke_params(); -pred = load_predicates(plant); % T_standby and FRET-predicate concretizations -figdir = fullfile('..', 'docs', 'figures'); -if ~exist(figdir, 'dir'), mkdir(figdir); end - -CtoF = @(T) T*9/5 + 32; - -%% ===== IC helpers ===== -% Operating IC: full-power steady state. -x0_op = pke_initial_conditions(plant); - -% Hot-standby shutdown IC: coolant flat at T_standby = T_c0 - 60 F ~ 275 C, -% reactor deep subcritical with only a trace neutron population. T_standby -% comes from reachability/predicates.json (single source of truth). Well -% inside the model's +/-50 C trust region around the operating point, and -% above coolant saturation at reduced operating pressure. -T_standby = pred.constants.T_standby; -n_shut = 1e-6; -C_shut = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_shut; -x0_shut = [n_shut; C_shut; T_standby; T_standby; T_standby]; - -% Heatup IC: reactor already taken critical at 0.1% power at hot-standby -% temperature. Mirrors real plant practice: achieve criticality, then -% heat up. Same T_standby as the shutdown IC — heatup begins from where -% shutdown left off. -n_heat = 1e-3; -C_heat = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_heat; -x0_heat = [n_heat; C_heat; T_standby; T_standby; T_standby]; - -%% ===== Mode 1: SHUTDOWN ===== -fprintf('\n===== Mode 1: ctrl_shutdown =====\n'); -Q_sg_shut = @(t) 0; % no SG demand -tspan_shut = [0, 600]; -[t1, X1, U1] = pke_solver(plant, Q_sg_shut, @ctrl_shutdown, [], tspan_shut, x0_shut); - -%% ===== Mode 2: HEATUP ===== -fprintf('\n===== Mode 2: ctrl_heatup =====\n'); -ref_heatup = struct(); -ref_heatup.T_start = T_standby; % ~275 C (hot standby, = T_c0 - 60 F) -ref_heatup.T_target = plant.T_c0; % 308.35 C (operating setpoint) -ref_heatup.ramp_rate = 28 / 3600; % 28 C/hr, tech-spec limit -Q_sg_heat = @(t) 0; % no SG demand during heatup -tspan_heat = [0, 5400]; % ~90 min (33 C span at 28 C/hr = 71 min + settling) -[t2, X2, U2] = pke_solver(plant, Q_sg_heat, @ctrl_heatup, ref_heatup, tspan_heat, x0_heat); - -%% ===== Mode 3a: OPERATION (plain P) ===== -fprintf('\n===== Mode 3a: ctrl_operation (P) =====\n'); -Q_sg_op = @(t) plant.P0 * (1.0 - 0.2 * (t >= 30)); % 100% -> 80% step -ref_op = struct('T_avg', plant.T_c0); -tspan_op = [0, 600]; -[t3, X3, U3] = pke_solver(plant, Q_sg_op, @ctrl_operation, ref_op, tspan_op, x0_op); - -%% ===== Mode 3b: OPERATION (LQR) ===== -fprintf('\n===== Mode 3b: ctrl_operation_lqr =====\n'); -[t3b, X3b, U3b] = pke_solver(plant, Q_sg_op, @ctrl_operation_lqr, [], tspan_op, x0_op); - -%% ===== Mode 4: SCRAM ===== -fprintf('\n===== Mode 4: ctrl_scram =====\n'); -% After a scram signal, the turbine trips and SG isolation occurs fast. -% Q_sg drops from P0 to a decay-heat-level sink (3% of P0) in ~10 s, -% then holds. Not a true decay-heat model — good enough to show the -% post-scram thermal trajectory without coolant going unphysically cold. -Q_sg_scr = @(t) plant.P0 * max(0.03, 1 - max(0, t - 10) / 10); -tspan_scr = [0, 600]; -[t4, X4, U4] = pke_solver(plant, Q_sg_scr, @ctrl_scram, [], tspan_scr, x0_op); - -%% ===== Per-mode 4-panel plots ===== -plot_pke_results(t1, X1, U1, plant, Q_sg_shut, 'ctrl\_shutdown (cold IC)'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_1_shutdown.png'), 'Resolution', 150); - -plot_pke_results(t2, X2, U2, plant, Q_sg_heat, 'ctrl\_heatup (ramp T\_avg)'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_2_heatup.png'), 'Resolution', 150); - -plot_pke_results(t3, X3, U3, plant, Q_sg_op, 'ctrl\_operation (P on T\_avg)'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3_operation.png'), 'Resolution', 150); - -plot_pke_results(t3b, X3b, U3b, plant, Q_sg_op, 'ctrl\_operation\_lqr'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3b_operation_lqr.png'), 'Resolution', 150); - -plot_pke_results(t4, X4, U4, plant, Q_sg_scr, 'ctrl\_scram'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_4_scram.png'), 'Resolution', 150); - -%% ===== Heatup ramp-tracking figure ===== -% Overlay the T_ref signal on T_avg so the lag is visible. -T_ref_trace = min(ref_heatup.T_start + ref_heatup.ramp_rate .* t2, ref_heatup.T_target); -figure('Position', [100 80 900 400], 'Name', 'Heatup ramp tracking'); -plot(t2/60, CtoF(T_ref_trace), 'k--', 'LineWidth', 1.2); hold on; -plot(t2/60, CtoF(X2(:,9)), 'r-', 'LineWidth', 1.5); -xlabel('Time [min]'); ylabel('T_{avg} [F]'); -legend('T_{ref}(t) (28 C/hr ramp)', 'T_{avg} (plant)', 'Location', 'southeast'); -title('ctrl\_heatup: reference tracking'); -grid on; -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_heatup_tracking.png'), 'Resolution', 150); - -%% ===== P vs LQR head-to-head on T_avg ===== -figure('Position', [100 80 900 400], 'Name', 'Operation: P vs LQR'); -plot(t3, CtoF(X3(:,9)), 'b-', 'LineWidth', 1.5); hold on; -plot(t3b, CtoF(X3b(:,9)), 'r-', 'LineWidth', 1.5); -yline(CtoF(plant.T_c0), 'k--'); -xlabel('Time [s]'); ylabel('T_{avg} [F]'); -legend('ctrl\_operation (P)', 'ctrl\_operation\_lqr', 'setpoint', 'Location', 'best'); -title('Operation mode: P vs LQR under 100% \rightarrow 80% Q_{sg} step'); -grid on; -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_op_P_vs_LQR.png'), 'Resolution', 150); - -%% ===== Summary figure: normalized power across all modes ===== -figure('Position', [100 80 1000 450], 'Name', 'Normalized power, all modes'); -subplot(2,2,1); semilogy(t1, max(X1(:,1), 1e-20)); grid on; - title('Shutdown'); xlabel('t [s]'); ylabel('n (log)'); -subplot(2,2,2); plot(t2/60, X2(:,1)); grid on; - title('Heatup'); xlabel('t [min]'); ylabel('n'); -subplot(2,2,3); plot(t3, X3(:,1)); grid on; - title('Operation'); xlabel('t [s]'); ylabel('n'); -subplot(2,2,4); semilogy(t4, max(X4(:,1), 1e-20)); grid on; - title('Scram'); xlabel('t [s]'); ylabel('n (log)'); -sgtitle('Normalized reactor power n(t) across DRC modes'); -exportgraphics(gcf, fullfile(figdir, 'mode_sweep_power_overview.png'), 'Resolution', 150); - -fprintf('\n=== Figures written to %s ===\n', figdir); diff --git a/plant-model/pke_initial_conditions.m b/plant-model/pke_initial_conditions.m deleted file mode 100644 index 74573fa..0000000 --- a/plant-model/pke_initial_conditions.m +++ /dev/null @@ -1,12 +0,0 @@ -function x0 = pke_initial_conditions(plant) -% PKE_INITIAL_CONDITIONS Steady-state initial condition vector. -% -% Returns x0 = [n0; C1..C6; T_f0; T_c0; T_cold0] (10 x 1) -% Precursor concentrations are the analytic equilibrium values at n = 1. - - n0 = 1.0; - C0 = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n0; - - x0 = [n0; C0; plant.T_f0; plant.T_c0; plant.T_cold0]; - -end diff --git a/plant-model/pke_linearize.m b/plant-model/pke_linearize.m deleted file mode 100644 index 7eb60fc..0000000 --- a/plant-model/pke_linearize.m +++ /dev/null @@ -1,54 +0,0 @@ -function [A, B, B_w, x_star, u_star, Q_star] = pke_linearize(plant, x_star, u_star, Q_star) -% PKE_LINEARIZE Numerical Jacobians of the PKE+T/H RHS at an operating point. -% -% Computes (A, B, B_w) such that, for small perturbations dx, du, dw -% around (x_star, u_star, Q_star): -% -% dx/dt ~ A*dx + B*du + B_w*dw, where w = Q_sg. -% -% Uses central finite differences. Step sizes are scaled to each -% variable's magnitude; fall back to an absolute epsilon near zero. -% -% Inputs: -% plant - parameter struct from pke_params() -% x_star - (optional) 10x1 linearization point; default = operating x0 -% u_star - (optional) scalar control; default = 0 -% Q_star - (optional) scalar SG heat removal [W]; default = plant.P0 -% -% Outputs: -% A - 10x10 state Jacobian df/dx -% B - 10x1 input Jacobian df/du -% B_w - 10x1 disturbance Jacobian df/d(Q_sg) -% x_star, u_star, Q_star - echoed back for the caller's records - - if nargin < 2 || isempty(x_star), x_star = pke_initial_conditions(plant); end - if nargin < 3 || isempty(u_star), u_star = 0; end - if nargin < 4 || isempty(Q_star), Q_star = plant.P0; end - - n = numel(x_star); - eps_rel = 1e-6; - eps_abs = 1e-8; - - % Wrap the RHS as a function of (x, u, w) with plant fixed. - f = @(x, u, w) pke_th_rhs(0, x, plant, @(t) w, u); - - f0 = f(x_star, u_star, Q_star); %#ok % sanity slot - - % --- A: df/dx via central differences, column by column --- - A = zeros(n, n); - for k = 1:n - h = max(eps_rel * abs(x_star(k)), eps_abs); - xp = x_star; xp(k) = xp(k) + h; - xm = x_star; xm(k) = xm(k) - h; - A(:, k) = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) / (2*h); - end - - % --- B: df/du --- - h = max(eps_rel * abs(u_star), eps_abs); - B = (f(x_star, u_star + h, Q_star) - f(x_star, u_star - h, Q_star)) / (2*h); - - % --- B_w: df/dQ_sg --- - h = max(eps_rel * abs(Q_star), 1.0); % Q_star ~ 1e9, abs floor at 1 W - B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) / (2*h); - -end diff --git a/plant-model/pke_params.m b/plant-model/pke_params.m deleted file mode 100644 index 1d6b85b..0000000 --- a/plant-model/pke_params.m +++ /dev/null @@ -1,86 +0,0 @@ -function p = pke_params() -% PKE_PARAMS Returns the parameter struct for the coupled PKE + T/H model. -% -% All plant parameters, kinetic data, and steady-state conditions live here. -% Modify this file to change the reactor you're modeling. - - p = struct(); - - % --- Neutronics --- - % Prompt neutron generation time for a typical PWR with enriched UO2 fuel - % and a water-moderated thermal spectrum. Literature range: 1e-5 to 5e-4 s; - % 1e-4 is a standard mid-range value (Duderstadt & Hamilton, Ch. 6). - p.Lambda = 1e-4; % prompt neutron generation time [s] - - % 6-group delayed neutron data for U-235 thermal fission (Keepin et al., - % Phys. Rev. 107, 1957). These are the canonical values used in virtually - % all PKE models. Total beta_eff ~ 650 pcm. - p.beta_i = [0.000215; 0.001424; 0.001274; 0.002568; 0.000748; 0.000273]; - p.lambda_i = [0.0124; 0.0305; 0.111; 0.301; 1.14; 3.01]; - p.beta = sum(p.beta_i); - - % --- Thermal-hydraulic --- - % 1000 MWth is a round number representative of a 2-loop PWR (e.g., - % Westinghouse 2-loop class). Actual plants range ~1800-3400 MWth for - % 2-loop to 4-loop designs. - p.P0 = 1000e6; % nominal thermal power [W] - - % Fuel mass: representative of total UO2 inventory in the core. - % A typical 4-loop PWR has ~100,000 kg UO2; scaled down for 1000 MWth. - p.M_f = 50000; % fuel mass [kg] - - % UO2 specific heat at operating temperatures (~400-800 C). Published - % values range 250-350 J/kg-K (MATPRO / Fink correlation); 300 is a - % standard lumped value for transient analysis. - p.c_f = 300; % fuel specific heat [J/kg-K] - - % Coolant mass in the core region. Based on core volume fraction of - % water in a PWR (~40% moderator-to-fuel ratio) and active core volume. - p.M_c = 20000; % coolant mass in core [kg] - - % Pressurized water specific heat at ~15.5 MPa, ~300 C. NIST/IAPWS - % gives cp ~ 5.4-5.5 kJ/kg-K in this range. - p.c_c = 5450; % coolant specific heat [J/kg-K] - - % Lumped fuel-to-coolant heat transfer coefficient x area. Chosen so that - % the steady-state fuel-to-coolant dT = P0/hA = 20 C, which is in the - % right ballpark for the average pellet-to-bulk-coolant temperature drop - % in a lumped single-node fuel model. - p.hA = 5e7; % fuel-to-coolant heat transfer coeff * area [W/K] - - % Primary coolant mass flow rate. Typical PWR: ~4000-6000 kg/s per loop. - % 5000 kg/s gives core dT = P0/(W*c_c) ~ 36.7 C, consistent with a - % ~35-38 C core rise in operating PWRs. - p.W = 5000; % coolant mass flow rate [kg/s] - - % Coolant mass in SG + hot/cold leg piping. Represents the ex-core primary - % inventory that participates in the loop transport delay. Larger than M_c - % because the SG tubes and piping hold significant volume. - p.M_sg = 30000; % coolant mass in SG + hot/cold leg piping [kg] - - % --- Reactivity feedback coefficients --- - % Doppler coefficient: fuel temperature feedback, predominantly from - % U-238 resonance broadening. Typical BOL range: -2 to -3 pcm/K - % (Todreas & Kazimi). -2.5 pcm/K is a standard mid-cycle value. - p.alpha_f = -2.5e-5; % fuel (Doppler) [dk/k per K] - - % Moderator temperature coefficient: captures water density change with - % temperature and its effect on neutron moderation/absorption. Typical - % range: -5 to -50 pcm/K depending on boron concentration and burnup. - % -10 pcm/K is moderate (mid-cycle, moderate boron). - p.alpha_c = -1.0e-4; % coolant (moderator) [dk/k per K] - - % --- Steady-state (derived from energy balance at P0, n=1) --- - % Q_sg = P0, all dT/dt = 0 - % T_hot - T_cold = P0 / (W*c_c) - % T_c = (T_hot + T_cold) / 2 - % T_f = T_c + P0 / hA - % Cold leg temperature of 290 C (554 F) is typical for a PWR at full - % power (Westinghouse plants: ~288-293 C depending on design). - p.T_cold0 = 290; % inlet coolant [C] - p.dT_core = p.P0 / (p.W * p.c_c); % core delta-T [C] - p.T_hot0 = p.T_cold0 + p.dT_core; - p.T_c0 = (p.T_hot0 + p.T_cold0) / 2; % avg coolant [C] - p.T_f0 = p.T_c0 + p.P0 / p.hA; % fuel [C] - -end diff --git a/plant-model/pke_solver.m b/plant-model/pke_solver.m deleted file mode 100644 index 12291d7..0000000 --- a/plant-model/pke_solver.m +++ /dev/null @@ -1,64 +0,0 @@ -function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan, x0) -% PKE_SOLVER Solve the coupled PKE + T/H system in closed loop with a mode controller. -% -% Inputs: -% plant - parameter struct from pke_params() -% Q_sg - function handle Q_sg(t) returning SG heat removal [W] -% ctrl_fn - function handle u = ctrl_fn(t, x, plant, ref) -% ref - struct of per-mode setpoints (e.g. ref.T_avg); [] if unused -% tspan - [t_start, t_end] in seconds -% x0 - (optional) initial state vector (10 x 1). Defaults to the -% operating steady state from pke_initial_conditions(plant). -% -% Outputs: -% t - time vector (M x 1) from the ODE solver -% X - state matrix (M x 10): columns are [n, C1..C6, T_f, T_c, T_cold] -% U - control input trajectory (M x 1): u evaluated at each (t(k), X(k,:)) - - CtoF = @(T) T * 9/5 + 32; - - %% ===== Print Steady-State ===== - fprintf('=== Steady-State Conditions ===\n'); - fprintf(' beta = %.5f (%.1f pcm)\n', plant.beta, plant.beta*1e5); - fprintf(' Lambda = %.1e s\n', plant.Lambda); - fprintf(' P0 = %.0f MWth\n', plant.P0/1e6); - fprintf(' T_cold = %.1f F\n', CtoF(plant.T_cold0)); - fprintf(' T_avg = %.1f F\n', CtoF(plant.T_c0)); - fprintf(' T_hot = %.1f F\n', CtoF(plant.T_hot0)); - fprintf(' T_fuel = %.1f F\n', CtoF(plant.T_f0)); - fprintf(' Core dT = %.1f F\n', plant.dT_core * 9/5); - fprintf('================================\n\n'); - - %% ===== Solve closed-loop ===== - if nargin < 6 || isempty(x0) - x0 = pke_initial_conditions(plant); - end - options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10); - - rhs = @(t, x) pke_th_rhs(t, x, plant, Q_sg, ctrl_fn(t, x, plant, ref)); - [t, X] = ode15s(rhs, tspan, x0, options); - - %% ===== Reconstruct control trajectory ===== - M = length(t); - U = zeros(M, 1); - for k = 1:M - U(k) = ctrl_fn(t(k), X(k,:).', plant, ref); - end - - %% ===== Print Final State ===== - n = X(end, 1); T_f = X(end, 8); T_c = X(end, 9); T_cold = X(end, 10); - T_hot = 2*T_c - T_cold; - rho_tot = U(end) ... - + plant.alpha_f*(T_f - plant.T_f0) ... - + plant.alpha_c*(T_c - plant.T_c0); - - fprintf('=== Final State (t = %.0f s) ===\n', t(end)); - fprintf(' Power = %.1f MWth (%.3f x nominal)\n', n*plant.P0/1e6, n); - fprintf(' T_fuel = %.1f F\n', CtoF(T_f)); - fprintf(' T_hot = %.1f F\n', CtoF(T_hot)); - fprintf(' T_avg = %.1f F\n', CtoF(T_c)); - fprintf(' T_cold = %.1f F\n', CtoF(T_cold)); - fprintf(' u = %.4f $ (external reactivity from controller)\n', U(end)/plant.beta); - fprintf(' rho = %.4f $ (total, incl. T-feedback)\n', rho_tot/plant.beta); - -end diff --git a/plant-model/pke_th_rhs.m b/plant-model/pke_th_rhs.m deleted file mode 100644 index f109224..0000000 --- a/plant-model/pke_th_rhs.m +++ /dev/null @@ -1,49 +0,0 @@ -function dxdt = pke_th_rhs(t, x, plant, Q_sg, u) -% PKE_TH_RHS ODE right-hand side for coupled PKE + thermal-hydraulics. -% -% State x = [n; C1..C6; T_f; T_c; T_cold] (10 x 1) -% -% Inputs: -% t - time [s] -% x - continuous state vector (10 x 1) -% plant - parameter struct from pke_params() -% Q_sg - function handle Q_sg(t) returning SG heat removal [W] -% u - external reactivity (rod worth) [dk/k], scalar -% -% The controller is responsible for computing u; this function treats u -% as a given input. Q_sg is the bounded disturbance. - - n = x(1); - C = x(2:7); - T_f = x(8); - T_c = x(9); - T_cold = x(10); - T_hot = 2 * T_c - T_cold; - - Q_sg_val = Q_sg(t); - - % --- Reactivity with feedback --- - % rho = u + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0) - rho = u ... - + plant.alpha_f * (T_f - plant.T_f0) ... - + plant.alpha_c * (T_c - plant.T_c0); - - % --- Neutronics --- - dndt = (rho - plant.beta) / plant.Lambda * n + plant.lambda_i' * C; - dCdt = (plant.beta_i / plant.Lambda) * n - plant.lambda_i .* C; - - % --- Thermal-hydraulics --- - % Fuel node: fission power in, heat transfer to coolant out - dTf_dt = (plant.P0 * n - plant.hA * (T_f - T_c)) / (plant.M_f * plant.c_f); - - % Core coolant: heat from fuel in, bulk enthalpy flow out - dTc_dt = (plant.hA * (T_f - T_c) - 2*plant.W*plant.c_c*(T_c - T_cold)) ... - / (plant.M_c * plant.c_c); - - % SG/loop node: hot leg enthalpy in, SG heat removal out - dTcold_dt = (plant.W * plant.c_c * (T_hot - T_cold) - Q_sg_val) ... - / (plant.M_sg * plant.c_c); - - dxdt = [dndt; dCdt; dTf_dt; dTc_dt; dTcold_dt]; - -end diff --git a/plant-model/plot_pke_results.m b/plant-model/plot_pke_results.m deleted file mode 100644 index 020712c..0000000 --- a/plant-model/plot_pke_results.m +++ /dev/null @@ -1,76 +0,0 @@ -function plot_pke_results(t, X, U, plant, Q_sg, plot_title) -% PLOT_PKE_RESULTS Standard 4-panel plot of PKE + T/H simulation results. -% -% Inputs: -% t - time vector from ODE solver -% X - state matrix (M x 10): [n, C1..C6, T_f, T_c, T_cold] -% U - control input trajectory (M x 1): scalar u at each time -% plant - parameter struct from pke_params() -% Q_sg - function handle Q_sg(t) used in the simulation -% plot_title - (optional) string used in the figure name and first subplot - - if nargin < 6 || isempty(plot_title) - plot_title = 'PKE + T/H closed-loop'; - end - - % --- Extract states --- - n = X(:, 1); - T_f = X(:, 8); - T_c = X(:, 9); - T_cold = X(:, 10); - T_hot = 2 * T_c - T_cold; - - % --- Convert temperatures to Fahrenheit for display --- - CtoF = @(T) T * 9/5 + 32; - T_f_F = CtoF(T_f); - T_c_F = CtoF(T_c); - T_cold_F = CtoF(T_cold); - T_hot_F = CtoF(T_hot); - - % --- Derived quantities --- - power_MW = n * plant.P0 / 1e6; - Q_sg_MW = arrayfun(Q_sg, t) / 1e6; - - rho_fb = plant.alpha_f * (T_f - plant.T_f0) + plant.alpha_c * (T_c - plant.T_c0); - rho_tot = U + rho_fb; - - % --- Plot --- - figure('Position', [100 50 1000 950], 'Name', plot_title); - - % Power & SG demand - subplot(4, 1, 1); - plot(t, power_MW, 'b-', 'LineWidth', 1.5); hold on; - plot(t, Q_sg_MW, 'r--', 'LineWidth', 1.2); - xlabel('Time [s]'); ylabel('Power [MW_{th}]'); - legend('Reactor Power', 'Q_{sg}', 'Location', 'east'); - title(sprintf('%s — Power and SG Heat Demand', plot_title)); - grid on; - - % Reactivity: controller u, feedback, total - subplot(4, 1, 2); - plot(t, U / plant.beta, 'k-', 'LineWidth', 1.5); hold on; - plot(t, rho_fb / plant.beta, 'b-.', 'LineWidth', 1.2); - plot(t, rho_tot / plant.beta, 'm:', 'LineWidth', 1.5); - xlabel('Time [s]'); ylabel('\rho [$]'); - legend('u (controller)', 'Feedback', 'Total', 'Location', 'east'); - title('Reactivity'); - grid on; - - % Coolant temperatures - subplot(4, 1, 3); - plot(t, T_hot_F, 'r-', 'LineWidth', 1.5); hold on; - plot(t, T_c_F, 'm--', 'LineWidth', 1.2); - plot(t, T_cold_F, 'b-', 'LineWidth', 1.5); - xlabel('Time [s]'); ylabel('Temperature [F]'); - legend('T_{hot}', 'T_{avg}', 'T_{cold}', 'Location', 'east'); - title('Coolant Temperatures'); - grid on; - - % Fuel temperature - subplot(4, 1, 4); - plot(t, T_f_F, 'Color', [0.85 0.33 0.1], 'LineWidth', 1.5); - xlabel('Time [s]'); ylabel('Temperature [F]'); - title('Fuel Temperature'); - grid on; - -end diff --git a/plant-model/test_heatup_rate.m b/plant-model/test_heatup_rate.m deleted file mode 100644 index 9978d72..0000000 --- a/plant-model/test_heatup_rate.m +++ /dev/null @@ -1,20 +0,0 @@ -%% test_heatup_rate.m — measure the actual heatup-rate trajectory from ctrl_heatup -clear; addpath('controllers'); addpath('../reachability'); -plant = pke_params(); -pred = load_predicates(plant); -x0 = [1e-3; (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * 1e-3; ... - pred.constants.T_standby; pred.constants.T_standby; pred.constants.T_standby]; -ref = struct('T_start', pred.constants.T_standby, ... - 'T_target', plant.T_c0, ... - 'ramp_rate', 28/3600); -[t, X, U] = pke_solver(plant, @(t) 0, @ctrl_heatup, ref, [0 5400], x0); - -a_f = 0.4587; a_c = -0.9587; a_cold = 0.5000; -dTcdt = a_f*X(:,8) + a_c*X(:,9) + a_cold*X(:,10); % C/s - -fprintf('\n=== Heatup rate statistics ===\n'); -fprintf(' max dT_c/dt: %.4f C/s = %6.1f C/hr\n', max(dTcdt), max(dTcdt)*3600); -fprintf(' min dT_c/dt: %.4f C/s = %6.1f C/hr\n', min(dTcdt), min(dTcdt)*3600); -fprintf(' rate limit: +/- 0.01389 C/s = +/- 50.0 C/hr\n'); -fprintf(' violates upper? %s\n', string(any(dTcdt > 0.01389))); -fprintf(' violates lower? %s\n', string(any(dTcdt < -0.01389))); diff --git a/plant-model/test_linearize.m b/plant-model/test_linearize.m deleted file mode 100644 index 3f99ffd..0000000 --- a/plant-model/test_linearize.m +++ /dev/null @@ -1,106 +0,0 @@ -%% test_linearize.m — sanity-check the Jacobians against the nonlinear sim -% -% Simulate two parallel trajectories under a small Q_sg step (5% down): -% 1. Full nonlinear plant with u = 0 -% 2. Linear dx/dt = A*dx + B_w*dw propagated by ode45 -% Compare deviations from x_op. For a small-enough disturbance, the -% trajectories should overlay to within a few percent over a few minutes. - -clear; clc; close all; -addpath('controllers'); - -plant = pke_params(); -x_op = pke_initial_conditions(plant); - -[A, B, B_w, xs, us, Qs] = pke_linearize(plant, x_op, 0, plant.P0); - -%% ===== Eigenvalue/conditioning check ===== -eigs_A = eig(A); -fprintf('\n=== Linearization at operating point ===\n'); -fprintf(' ||A|| = %.3e\n', norm(A)); -fprintf(' ||B|| = %.3e\n', norm(B)); -fprintf(' ||B_w|| = %.3e\n', norm(B_w)); -fprintf(' stable? = %s (max Re(eig) = %.3e)\n', ... - iif(all(real(eigs_A) < 1e-8), 'YES', 'NO'), max(real(eigs_A))); -fprintf(' fastest pole ~ %.3e /s\n', max(abs(real(eigs_A)))); -fprintf(' slowest pole ~ %.3e /s\n', min(abs(real(eigs_A(real(eigs_A) < -1e-10))))); - -%% ===== Linear vs nonlinear trajectory under small Q_sg step ===== -% 5% down step at t = 30s. -dQ = -0.05 * plant.P0; -Q_sg = @(t) plant.P0 + dQ * (t >= 30); - -tspan = [0, 300]; - -% Nonlinear -opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10); -rhs_nl = @(t, x) pke_th_rhs(t, x, plant, Q_sg, 0); -[t_nl, X_nl] = ode15s(rhs_nl, tspan, x_op, opts); - -% Linear (deviation) -dQ_of_t = @(t) dQ * (t >= 30); -rhs_lin = @(t, dx) A*dx + B_w*dQ_of_t(t); -[t_lin, DX_lin] = ode15s(rhs_lin, tspan, zeros(size(x_op)), opts); - -% Reconstruct absolute linear state for plotting -X_lin = DX_lin + x_op.'; - -%% ===== Plot comparison — focus on thermal states ===== -CtoF = @(T) T*9/5 + 32; -figdir = fullfile('..', 'docs', 'figures'); -if ~exist(figdir, 'dir'), mkdir(figdir); end - -figure('Position', [100 80 1100 700], 'Name', 'Linear vs nonlinear sanity'); - -subplot(2,2,1); -plot(t_nl, X_nl(:,1), 'b-', 'LineWidth', 1.5); hold on; -plot(t_lin, X_lin(:,1), 'r--','LineWidth', 1.5); -xlabel('t [s]'); ylabel('n'); grid on; -title('Normalized power'); legend('nonlinear', 'linear', 'Location', 'best'); - -subplot(2,2,2); -plot(t_nl, CtoF(X_nl(:,8)), 'b-', 'LineWidth', 1.5); hold on; -plot(t_lin, CtoF(X_lin(:,8)), 'r--','LineWidth', 1.5); -xlabel('t [s]'); ylabel('T_f [F]'); grid on; -title('Fuel temperature'); - -subplot(2,2,3); -plot(t_nl, CtoF(X_nl(:,9)), 'b-', 'LineWidth', 1.5); hold on; -plot(t_lin, CtoF(X_lin(:,9)), 'r--','LineWidth', 1.5); -xlabel('t [s]'); ylabel('T_c [F]'); grid on; -title('Avg coolant'); - -subplot(2,2,4); -plot(t_nl, CtoF(X_nl(:,10)), 'b-', 'LineWidth', 1.5); hold on; -plot(t_lin, CtoF(X_lin(:,10)),'r--','LineWidth', 1.5); -xlabel('t [s]'); ylabel('T_{cold} [F]'); grid on; -title('Cold-leg coolant'); - -sgtitle('Linear-model sanity check: 5% Q_{sg} down-step, u = 0'); -exportgraphics(gcf, fullfile(figdir, 'linearize_sanity.png'), 'Resolution', 150); - -%% ===== Quantitative error at a few times ===== -fprintf('\n=== Max |linear - nonlinear| at sampled times ===\n'); -fprintf(' (normalized by operating-point magnitude where nonzero)\n'); -for ts = [60, 120, 300] - xi_nl = interp1(t_nl, X_nl, ts); - xi_lin = interp1(t_lin, X_lin, ts); - rel_err = abs(xi_nl - xi_lin) ./ max(abs(x_op.'), 1e-6); - fprintf(' t = %3d s: max rel err = %.2e (%s)\n', ts, max(rel_err), ... - state_name(find(rel_err == max(rel_err), 1))); -end - -%% ===== Save linearization ===== -save(fullfile('..', 'reachability', 'linearization_at_op.mat'), ... - 'A', 'B', 'B_w', 'xs', 'us', 'Qs', '-v7'); -fprintf('\nSaved A/B/B_w to ../reachability/linearization_at_op.mat\n'); - -%% ===== helpers ===== -function s = state_name(k) - names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'}; - s = names{k}; -end - -function y = iif(cond, a, b) - if cond, y = a; else, y = b; end -end diff --git a/reachability/README.md b/reachability/README.md index 520391f..10bb4cc 100644 --- a/reachability/README.md +++ b/reachability/README.md @@ -4,95 +4,95 @@ Continuous-mode verification for the PWR_HYBRID_3 hybrid controller. ## Soundness status: APPROXIMATE -The current `reach_operation.m` result is **not a sound reach tube for -the physical plant**. It is a sound over-approximation of the -*linearized* closed-loop system (A_cl = A - BK around x_op) under -bounded disturbance. The linear model is itself an approximation of -the nonlinear plant (`../plant-model/pke_th_rhs.m`), and that -approximation error is not currently bounded or inflated into the tube. +The current linear-reach results are **not a sound reach tube for +the physical plant**. They are sound over-approximations of the +*linearized* closed-loop system ($A_{\mathrm{cl}} = A - BK$ around +$x_{\mathrm{op}}$) under bounded disturbance. The linear model is +itself an approximation of the nonlinear plant (`../code/src/pke_th_rhs.jl`), +and that approximation error is not currently bounded or inflated +into the tube. Two paths to upgrade to a sound result: -1. **Nonlinear reach directly** — CORA `nonlinearSys`, JuliaReach - `BlackBoxContinuousSystem`, or equivalent. More expensive but the - honest answer. +1. **Nonlinear reach directly** — TMJets (Taylor-model integration) + via ReachabilityAnalysis.jl. Currently limited to ~10-second + horizons by prompt-neutron stiffness; needs a reduced-order PKE + (prompt-jump approximation) to extend to mode-obligation horizons. 2. **Linear reach + Taylor-remainder inflation** — compute an upper - bound on `||f_nl(x, u) - (A x + B u)||` over the reach set (via - Hessian norm estimate on each component of `f_nl`) and inflate the - linear tube by that bound. Less expensive, still rigorous. + bound on `||f_nl(x, u) - (A x + B u)||` over the reach set and + inflate the linear tube by that bound. Cheaper, still rigorous. -Both are thesis-blocking for any safety claim. Deferred only until -the per-mode plumbing is solid; it is not a "nice to have". - -The current 5-orders-of-margin buffer (reach envelope ~0.03 K against -a 5 K safety band) means linearization error would have to be huge to -invalidate the conclusion, but that is vibes, not a proof. +Both are thesis-blocking for any safety claim. The current +5-orders-of-margin buffer (reach envelope ~0.03 K against a +5 K safety band on $T_c$) means linearization error would have to +be huge to invalidate the conclusion — but that's vibes, not a proof. ## Related open issues -- **Saturation semantics.** `ctrl_heatup.m` uses `sat(u, u_min, u_max)`. - Saturation is formally a 3-mode piecewise-affine system. For - heatup reach this has to be handled as (a) hybrid locations, or - (b) proven dormant via reach on `u_unsat`. Not modeled in the +- **Saturation semantics.** `ctrl_heatup.jl` uses `clamp(u, u_min, u_max)`. + Saturation is formally a 3-mode piecewise-affine system. For + heatup reach this must be handled as (a) hybrid locations, or + (b) proven dormant via reach on `u_unsat`. Not modeled in the current artifacts (operation-mode LQR saturation is dormant in practice but the proof is implicit). -- **Parametric uncertainty in α_f, α_c.** Real plants have α drift - with burnup (~20%), boron (α_c ranges 10×), xenon. The - feedback-linearization in `ctrl_heatup.m` assumes exact α; a robust +- **Parametric uncertainty in α_f, α_c.** Real plants have α drift + with burnup (~20%), boron (α_c ranges 10×), xenon. The + feedback-linearization in `ctrl_heatup` assumes exact α; a robust treatment would make α an interval and propagate parametric reach. - Currently idealized — flag in the chapter. + Currently idealized. ## What's here -**Per-mode only.** Following the compositionality argument in the thesis: -verify each continuous mode separately, let the DRC handle discrete -switching. Current focus: **operation mode** under LQR feedback. +**Per-mode only.** Following the compositionality argument in the +thesis: verify each continuous mode separately, let the DRC handle +discrete switching. -## What's here +Files: -- `linearization_at_op.mat` — A, B, B_w and reference point, generated by - `../plant-model/test_linearize.m`. -- `reach_linear.m` — box-zonotope propagation of the closed-loop linear - model under bounded disturbance. Pure MATLAB, no external toolbox. -- `barrier_lyapunov.m` — Lyapunov-ellipsoid barrier certificate for the - closed-loop linear system. Solves a Lyapunov equation, reports the - smallest sub-level set containing the initial set and closed under - the disturbance. -- `reach_operation.m` — end-to-end operation-mode reach: linearize at - x_op, compute LQR gain, propagate zonotope reach set, check against - the `t_avg_in_range` predicate. -- `figures/` — generated plots. +- `predicates.json` — single source of truth for predicate + concretizations. Three groups: + - `operational_deadbands` — soft bands used by DRC for mode + transitions (`t_avg_above_min`, `t_avg_in_range`, `p_above_crit`). + - `safety_limits` — hard one-sided halfspaces (fuel centerline, + trip setpoints, subcooling, heatup-rate bounds). + - `mode_invariants` — `inv1_holds`, `inv2_holds` as conjunctions + of safety limits. + - `mode_boundaries` — per-mode $X_{\mathrm{entry}}$, $X_{\mathrm{safe}}$, + $X_{\mathrm{exit}}$, $T_{\min}$, $T_{\max}$. + +- `WALKTHROUGH.md` — standalone document explaining the reach-obligation + taxonomy, per-mode entry/exit definitions, current results, soundness + status. Read this for the full story. + +The reach code itself lives in `../code/`: + +- `../code/src/reach_linear.jl` — hand-rolled box reach propagator. +- `../code/src/load_predicates.jl` — reads `predicates.json`. +- `../code/scripts/reach_operation.jl` — operation-mode linear reach. +- `../code/scripts/barrier_lyapunov.jl` — Lyapunov barrier attempt. +- `../code/scripts/barrier_compare_OL_CL.jl` — OL vs CL comparison. +- `../code/scripts/reach_heatup_nonlinear.jl` — TMJets nonlinear reach + (heatup, saturation-disabled, 10-second horizon cap). ## Running -From MATLAB: - -```matlab -cd reachability -reach_operation % computes reach set + plots -barrier_lyapunov % solves Lyapunov, reports invariant ellipsoid +```bash +cd code +julia --project=. scripts/reach_operation.jl +julia --project=. scripts/barrier_lyapunov.jl +julia --project=. scripts/barrier_compare_OL_CL.jl +julia --project=. scripts/reach_heatup_nonlinear.jl ``` -## Tool choice - -Currently using a hand-rolled zonotope reach because: -- Avoids a ~0.5 GB CORA install for a first-pass result. -- Linear reach with bounded disturbance has a clean analytic form - (matrix exponential on the state, integral of e^(A(t-s))·B_w·w ds - for the disturbance). -- Stays inside MATLAB, which is where the plant model lives. - -If we need nonlinear reach (and we will, for non-LQR controllers or -larger reach sets where linearization error matters), the planned -options are CORA (MATLAB) or JuliaReach (port the plant to Julia). +Results save here as `*.mat` (gitignored). ## What this does NOT do yet - Any sound reach tube (see top of this file). -- Nonlinear reach for the original P controller on operation. -- Heatup reach (ramped reference makes x* time-varying — needs - trajectory-LQR or a different formulation, and the saturation - semantics need to be made explicit). +- Nonlinear reach on horizons > 10 s (needs reduced-order PKE). - Shutdown, scram, initialization reach. - Hybrid-system level verification (mode switching validity). - Parametric robustness to α_f, α_c drift. +- Polytopic or SOS barriers — the canonical quadratic Lyapunov + barrier fails structurally on this plant (see `WALKTHROUGH.md` + and the OL-vs-CL comparison script). diff --git a/reachability/WALKTHROUGH.md b/reachability/WALKTHROUGH.md index b8a8b66..05a74f6 100644 --- a/reachability/WALKTHROUGH.md +++ b/reachability/WALKTHROUGH.md @@ -128,7 +128,13 @@ will require real tech-spec numbers. ## 3. Code walkthrough -Files in `reachability/` and what each does. +Files in `code/` (Julia) and `reachability/` (JSON + docs) and what each does. + +**Note:** as of 2026-04-20, all reach / plant / controller code is Julia. +The pipeline was ported end-to-end from an earlier MATLAB implementation. +Numerical results match MATLAB to machine precision on per-halfspace +margins; trajectories match to ~3 decimals (different ODE solvers, +same stiff family). ### `predicates.json` @@ -146,7 +152,7 @@ categories: Plus `mode_boundaries` (entry/exit/time budgets per mode, as above). -### `load_predicates.m` +### `code/src/load_predicates.jl` Reads the JSON, resolves `rhs_expr` strings (which reference plant-derived constants like `T_c0`, `T_cold0`, `T_standby`) into numeric @@ -163,7 +169,7 @@ pred.mode_invariants.inv2_holds % conjunction of safety_limits Each `A_poly` is an `n_halfspaces x 10` matrix; each `b_poly` is `n_halfspaces x 1`. The predicate is `{x : A_poly * x ≤ b_poly}`. -### `pke_linearize.m` (in `plant-model/`) +### `code/src/pke_linearize.jl` Numerical Jacobians via central finite differences at any `(x_star, u_star, Q_star)`. Returns `(A, B, B_w)` such that for small @@ -173,7 +179,7 @@ perturbations `dx`, `du`, `dw`: dx/dt ≈ A dx + B du + B_w dw, w = Q_sg ``` -Validated in `plant-model/test_linearize.m`: under a 5% Q_sg step at +Validated in `code/scripts/test_linearize.jl` (via `sim_sanity.jl`): under a 5% Q_sg step at `u=0`, the linear model matches the nonlinear sim to ~4e-4 relative error at 60 s, improving to 5e-6 at 300 s. @@ -189,7 +195,7 @@ See `docs/figures/linearize_sanity.png`: ![Linear vs nonlinear sanity](../docs/figures/linearize_sanity.png) -### `reach_linear.m` +### `code/src/reach_linear.jl` ~50 lines. Propagates a box over-approximation of the reach tube for a linear system `dx/dt = A_cl·x + B_w·w`, with `x(0)` in a box and `w` in @@ -225,7 +231,7 @@ M_aug = expm([A_cl, B_w; zeros(1, n+1)] * dt); G_step = M_aug(1:n, n+1); % upper-right block IS the integral ``` -### `reach_operation.m` +### `code/scripts/reach_operation.jl` End-to-end operation-mode reach. Pipeline: @@ -255,9 +261,9 @@ for each halfspace a' x <= b in inv2_holds: Reports per-halfspace margins so failure modes are attributable. -### `barrier_lyapunov.m` +### `code/scripts/barrier_lyapunov.jl` -Analytic counterpart to `reach_operation.m`. Solves a Lyapunov equation +Analytic counterpart to `code/scripts/reach_operation.jl`. Solves a Lyapunov equation to get an invariant ellipsoid, checks whether the ellipsoid fits inside the safety limits. @@ -277,10 +283,10 @@ max |a' dx| on {dx : dx' P dx ≤ γ} = sqrt(γ · a' P^(-1) a) ``` Sweeps `Qbar(9,9)` from 10 to 1e6 to find the best quadratic barrier. -See `barrier_compare_OL_CL.m` for the open-loop vs closed-loop +See `code/scripts/barrier_compare_OL_CL.jl` for the open-loop vs closed-loop comparison. -### `barrier_compare_OL_CL.m` +### `code/scripts/barrier_compare_OL_CL.jl` Side-by-side: run the Lyapunov barrier on open-loop `A` (no controller) vs closed-loop `A_cl = A - BK`. Written to answer the question "is @@ -361,14 +367,14 @@ it's a structural issue not a tuning one. ### Heatup mode — simulation only, no reach yet -The `main_mode_sweep.m` heatup scenario drives the plant from hot +The `code/scripts/main_mode_sweep.jl` heatup scenario drives the plant from hot standby through the operating temperature: ![Heatup tracking](../docs/figures/mode_sweep_heatup_tracking.png) Controller tracks the 28 °C/hr ramp with ~2 °F lag and ~5 °F final overshoot. **Peak rate during ramp is 48.5 °C/hr** (verified by -`plant-model/test_heatup_rate.m`), right at the 50 °C/hr placeholder +`code/scripts/test_heatup_rate.jl`), right at the 50 °C/hr placeholder invariant but 70% over the nominal tech-spec rate. A strict 28 °C/hr invariant would be violated by the current controller tuning. @@ -405,7 +411,7 @@ Power overview: **The current reach tube is approximate, not sound, for the nonlinear plant.** -`reach_linear.m` produces a sound over-approximation of the *linearized* +`code/src/reach_linear.jl` produces a sound over-approximation of the *linearized* closed-loop system's reach set. The linear model approximates the nonlinear plant with an error that is bounded but not currently computed or propagated into the tube. Two paths to soundness: @@ -449,7 +455,7 @@ vibes, not a proof. ### What's next 1. **Nonlinear reach on heatup via JuliaReach — partial result in.** - `../julia-port/scripts/reach_heatup_nonlinear.jl` ran successfully + `code/scripts/reach_heatup_nonlinear.jl` ran successfully to T=10 s on the full 10-state nonlinear closed-loop (saturation disabled). Got 10,583 reach-sets with envelope T_c ∈ [274.45, 295] and n ∈ [-5e-4, 5e-3]. Nonlinear reach is *functional* in Julia — diff --git a/reachability/barrier_compare_OL_CL.m b/reachability/barrier_compare_OL_CL.m deleted file mode 100644 index 8147902..0000000 --- a/reachability/barrier_compare_OL_CL.m +++ /dev/null @@ -1,86 +0,0 @@ -%% barrier_compare_OL_CL.m — does LQR help the Lyapunov barrier, or not? -% -% Head-to-head: solve the same Lyapunov-ellipsoid barrier analysis on -% (i) the open-loop plant A (u = 0) -% (ii) the LQR closed loop A_cl = A - B*K -% with identical Qbar and identical disturbance bound. Report per- -% halfspace bounds against inv2_holds so we can see exactly where LQR -% helps and where it doesn't. -% -% The plant is already Hurwitz open-loop (max Re(eig) = -0.012 /s), so -% the Lyapunov equation has a solution in both cases. - -clear; clc; -addpath('../plant-model', '../plant-model/controllers'); - -plant = pke_params(); -x_op = pke_initial_conditions(plant); -pred = load_predicates(plant); - -[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); - -%% ===== LQR gain ===== -Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1]); -R_lqr = 1e6; -try - K = lqr(A, B, Q_lqr, R_lqr); -catch - [~, ~, K] = icare(A, B, Q_lqr, R_lqr); -end -A_cl = A - B*K; - -%% ===== Shared pieces ===== -Qbar = diag([1, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1, 1e2, 1]); % best from the earlier sweep -w_bar = 0.15 * plant.P0; - -delta_entry = [0.01 * x_op(1); - 0.001 * abs(x_op(2:7)); - 0.1; 0.1; 0.1]; - -inv2 = pred.mode_invariants.inv2_holds; -A_inv = inv2.A_poly; b_inv = inv2.b_poly; comps = inv2.components; -b_dev = b_inv - A_inv * x_op; % headroom in deviation coords - -%% ===== Helper ===== - function [gamma, maxima, eig_cl] = run_case(Acase, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev) - % Solve Lyapunov, compute invariant level, per-halfspace maxima. - P = lyap(Acase.', Qbar); - Ph = sqrtm(P); Phi = inv(Ph); - mu = min(eig(Phi * Qbar * Phi)); - g_bound = sqrt(B_w.' * P * B_w); - c_inv = (2 * w_bar * g_bound / mu)^2; - c_entry = delta_entry.' * abs(P) * delta_entry; - gamma = max(c_inv, c_entry); - Pinv = inv(P); - maxima = zeros(size(A_inv, 1), 1); - for k = 1:size(A_inv, 1) - a = A_inv(k, :).'; - maxima(k) = sqrt(gamma * (a.' * Pinv * a)); - end - eig_cl = eig(Acase); - end - -%% ===== Run both ===== -[gamma_OL, max_OL, eig_OL] = run_case(A, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev); -[gamma_CL, max_CL, eig_CL] = run_case(A_cl, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev); - -%% ===== Report ===== -fprintf('\n=== Open-loop vs LQR closed-loop Lyapunov barrier ===\n'); -fprintf(' max Re(eig) OL = %.3e CL = %.3e\n', max(real(eig_OL)), max(real(eig_CL))); -fprintf(' min Re(eig) OL = %.3e CL = %.3e\n', min(real(eig_OL)), min(real(eig_CL))); -fprintf(' gamma OL = %.3e CL = %.3e (ratio CL/OL = %.3g)\n', ... - gamma_OL, gamma_CL, gamma_CL/gamma_OL); - -fprintf('\n Per-halfspace max |a'' dx| on gamma-ellipsoid:\n'); -fprintf(' %-22s %-10s %-12s %-12s %-12s %-10s\n', ... - 'halfspace', 'headroom', 'open-loop', 'closed-loop', 'CL - OL', 'ratio'); -for k = 1:size(A_inv, 1) - ratio = max_CL(k) / max_OL(k); - fprintf(' %-22s %10.3f %12.3f %12.3f %+12.3f %8.3fx\n', ... - comps{k}, b_dev(k), max_OL(k), max_CL(k), max_CL(k) - max_OL(k), ratio); -end - -fprintf('\n Interpretation:\n'); -fprintf(' - CL < OL on a row ==> LQR tightens that halfspace.\n'); -fprintf(' - CL ~= OL ==> LQR does not help that direction at all.\n'); -fprintf(' - CL > OL ==> LQR made that direction WORSE for the barrier (!!).\n'); diff --git a/reachability/barrier_lyapunov.m b/reachability/barrier_lyapunov.m deleted file mode 100644 index da8a6d1..0000000 --- a/reachability/barrier_lyapunov.m +++ /dev/null @@ -1,200 +0,0 @@ -%% barrier_lyapunov.m — Lyapunov-ellipsoid barrier certificate -% -% For dx/dt = A_cl x + B_w w with A_cl Hurwitz and ||w||_inf <= w_bar: -% -% 1. Solve A_cl' P + P A_cl = -Qbar (Qbar > 0, chosen = I). -% Then V(x) = x' P x is a Lyapunov function for the undisturbed -% system, with dV/dt = -x'x - x'(Qbar-I)x (here Qbar=I gives -x'x). -% -% 2. Under bounded disturbance: -% dV/dt = -x'x + 2 x' P B_w w -% <= -||x||^2 + 2 ||P B_w|| w_bar ||x||. -% dV/dt <= 0 whenever ||x|| >= 2 ||P B_w|| w_bar. -% So the ball B_r := {x : ||x|| <= 2 ||P B_w|| w_bar} contains -% the set where V can still grow. Any level set {V <= c} that -% contains B_r is forward-invariant. -% -% 3. Smallest such c: c* = lambda_max(P) * r^2, where r = 2||P B_w||w_bar. -% -% 4. Safety: the barrier is B(x) = V(x) - gamma, with gamma chosen -% large enough to contain X_entry but small enough that the level -% set stays inside X_safe. We report whether such a gamma exists. -% -% This is an ellipsoidal over-approximation, generally much looser than -% the box/zonotope reach in reach_operation.m, but it gives a *certificate* -% (a closed-form invariant function) rather than just a numerical tube. - -clear; clc; -addpath('../plant-model', '../plant-model/controllers'); - -plant = pke_params(); -x_op = pke_initial_conditions(plant); - -%% ===== Build A_cl, B_w ===== -[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); - -Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e-2, 1e2, 1]); -R_lqr = 1e6; -try - K = lqr(A, B, Q_lqr, R_lqr); -catch - [~, ~, K] = icare(A, B, Q_lqr, R_lqr); -end -A_cl = A - B*K; - -%% ===== Solve Lyapunov equation ===== -% A_cl' P + P A_cl + Qbar = 0. Qbar shaped to weight T_c heavily so the -% resulting ellipsoidal invariant sets are tight in the T_c direction. -% Without shaping, isotropic Qbar = I gives ellipsoids stretched along -% the slow-precursor directions, making the T_c safety bound useless. -Qbar = diag([1, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1, 1e3, 1]); -P = lyap(A_cl.', Qbar); -assert(all(eig(P) > 0), 'P not positive definite'); - -%% ===== Safety spec (used by sweep and final check) ===== -% Load inv2_holds (conjunction of safety halfspaces) from the predicates -% source of truth. Each row k of A_inv defines a halfspace a_k' x <= b_k; -% the barrier must bound max(a_k' * dx) over the ellipsoid for each k. -addpath('../plant-model'); -pred = load_predicates(plant); -inv2 = pred.mode_invariants.inv2_holds; -A_inv = inv2.A_poly; -b_inv = inv2.b_poly; -comp_names = inv2.components; -n_halfspaces = size(A_inv, 1); - -% Convert limits to deviation from x_op: -% for halfspace a' x <= b, the deviation-frame bound is a' dx <= b - a' x_op. -b_inv_dev = b_inv - A_inv * x_op; - -% Backward-compat scalars for existing prints. -e9 = zeros(10, 1); e9(9) = 1; -delta_safe_Tc = 5.0; - -%% ===== Disturbance bound ===== -% |w| <= w_bar where w = Q_sg - Q_nom. Take the same 15% down-load as -% reach_operation.m. -w_bar = 0.15 * plant.P0; - -% --- Invariant-level computation --- -% dV/dt = -x' Qbar x + 2 x' P B_w w. -% Taking the worst w = w_bar * sign(x' P B_w), the scalar g = x' P B_w: -% dV/dt <= -x' Qbar x + 2 w_bar |g|. -% Let u = P^{1/2} x (so V = ||u||^2). Then |g| = |u' P^{-1/2} P B_w| -% <= ||u|| * ||P^{-1/2} P B_w|| = sqrt(V) * sqrt(B_w' P B_w). -% And x' Qbar x >= lambda_min(P^{-1/2} Qbar P^{-1/2}) * V (call this mu). -% So dV/dt <= -mu V + 2 w_bar sqrt(B_w' P B_w) sqrt(V). -% dV/dt <= 0 whenever sqrt(V) >= 2 w_bar sqrt(B_w' P B_w) / mu, -% i.e. V >= (2 w_bar sqrt(B_w' P B_w) / mu)^2 := c_inv. -% -% This is much tighter than the isotropic ball bound — it uses the fact -% that B_w only pokes one direction of the ellipsoid. -P_half = sqrtm(P); -P_half_inv = inv(P_half); -mu = min(eig(P_half_inv * Qbar * P_half_inv)); -g_bound = sqrt(B_w.' * P * B_w); % sqrt(B_w' P B_w) -c_inv = (2 * w_bar * g_bound / mu)^2; - -fprintf('\n=== Lyapunov barrier certificate ===\n'); -fprintf(' lambda_min(P) = %.3e\n', min(eig(P))); -fprintf(' lambda_max(P) = %.3e\n', max(eig(P))); -fprintf(' sqrt(B_w'' P B_w) = %.3e\n', g_bound); -fprintf(' mu (Qbar eig on P-metric) = %.3e\n', mu); -fprintf(' w_bar (15%% P0) = %.3e W\n', w_bar); -fprintf(' c_inv (invariant level) = %.3e\n', c_inv); - -%% ===== Containment of initial set ===== -% Initial set: box around x_op with halfwidth delta_entry (matches reach_operation). -delta_entry = [0.01 * x_op(1); - 0.001 * abs(x_op(2:7)); - 0.1; 0.1; 0.1]; - -% Worst-case V over the initial box: max x'Px = sum over all 2^n corners. -% For small n we could enumerate, but the sharper bound is -% max V(dx) = delta_entry' * |P| * delta_entry -% (elementwise abs of P), which is the L1 energy bound. -c_entry = delta_entry.' * abs(P) * delta_entry; - -fprintf('\n c_entry (bound on V over initial box) = %.3e\n', c_entry); - -gamma = max(c_entry, c_inv); % barrier level must contain both -fprintf(' gamma (barrier level) = %.3e\n', gamma); -if gamma == c_entry - fprintf(' (initial set drives gamma — invariant piece already inside entry)\n'); -else - fprintf(' (disturbance inflation drives gamma)\n'); -end - -%% ===== Sweep Qbar(9,9) to find the tightest safe barrier ===== -% The isotropic Lyapunov bound is conservative because the "slow decay" -% direction dominates mu even when T_c is tightly controlled. Sweep the -% T_c weight to find a Qbar that yields a sub-5K barrier if one exists -% for this LQR design. -fprintf('\n=== Sweeping Qbar(T_c) weight ===\n'); -weights = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6]; -best_dTc = inf; best_w = NaN; best_gamma = NaN; best_P = []; -for wTc = weights - Qbar_s = Qbar; Qbar_s(9,9) = wTc; - try - Ps = lyap(A_cl.', Qbar_s); - catch - continue - end - if any(eig(Ps) <= 0), continue, end - Ph = sqrtm(Ps); Phi = inv(Ph); - mu_s = min(eig(Phi * Qbar_s * Phi)); - g_s = sqrt(B_w.' * Ps * B_w); - ci_s = (2 * w_bar * g_s / mu_s)^2; - ce_s = delta_entry.' * abs(Ps) * delta_entry; - g_s_level = max(ci_s, ce_s); - Pinv_s = inv(Ps); - dTc_s = sqrt(g_s_level * (e9.' * Pinv_s * e9)); - fprintf(' Qbar(9,9) = %.0e -> gamma = %.3e, max|dT_c| = %7.3f K\n', ... - wTc, g_s_level, dTc_s); - if dTc_s < best_dTc - best_dTc = dTc_s; best_w = wTc; best_gamma = g_s_level; best_P = Ps; - end -end -fprintf(' Best: Qbar(9,9) = %.0e -> max|dT_c| = %.3f K\n', best_w, best_dTc); -if best_dTc <= delta_safe_Tc - fprintf(' *** TIGHT BARRIER FOUND: V(x) = x.'' P_best x - gamma ***\n'); - P = best_P; gamma = best_gamma; -end - -%% ===== Safety: does the gamma-level set fit inside X_safe? ===== -% X_safe = { x : |T_c - T_c0| <= 5 K }, i.e. |e_9.' * dx| <= 5. -% Max |e_9.' * dx| over {dx : dx' P dx <= gamma} is sqrt(gamma * e_9' P^-1 e_9). -Pinv = inv(P); -max_dTc_on_ellipsoid = sqrt(gamma * (e9.' * Pinv * e9)); - -fprintf('\n=== Safety check on T_c ===\n'); -fprintf(' Max |dT_c| on gamma-ellipsoid = %.3f K\n', max_dTc_on_ellipsoid); -fprintf(' Safe band = +/- %.1f K\n', delta_safe_Tc); -if max_dTc_on_ellipsoid <= delta_safe_Tc - fprintf(' BARRIER VALID: V(x) = x.''Px - %.3e certifies T_c safety.\n', gamma); -else - fprintf(' *** BARRIER TOO LOOSE *** - ellipsoid reach into unsafe region.\n'); - fprintf(' Try a tighter LQR (bigger Q_Tc or smaller R) or tighter X_entry.\n'); -end - -%% ===== Per-halfspace barrier check against inv2_holds ===== -% For each safety halfspace a' dx <= b_dev, the max of a' dx over the -% gamma-ellipsoid {dx : dx' P dx <= gamma} is sqrt(gamma * a' P^{-1} a). -% Compare to b_dev (the headroom from x_op to the safety limit). -fprintf('\n=== Lyapunov barrier vs inv2_holds halfspaces ===\n'); -Pinv = inv(P); -for k = 1:n_halfspaces - a = A_inv(k, :).'; - max_adx = sqrt(gamma * (a.' * Pinv * a)); - margin = b_inv_dev(k) - max_adx; - status = 'OK'; - if margin < 0, status = '*** BARRIER TOO LOOSE ***'; end - fprintf(' [%-20s] headroom = %8.3f | max a''dx = %8.3f | margin = %+8.3f %s\n', ... - comp_names{k}, b_inv_dev(k), max_adx, margin, status); -end - -save(fullfile('.', 'barrier_lyapunov_result.mat'), ... - 'P', 'gamma', 'c_entry', 'c_inv', 'w_bar', 'K', 'A_cl', 'delta_entry', ... - 'max_dTc_on_ellipsoid', 'delta_safe_Tc', 'A_inv', 'b_inv_dev', '-v7'); - -fprintf('\nSaved barrier to ./barrier_lyapunov_result.mat\n'); diff --git a/reachability/load_predicates.m b/reachability/load_predicates.m deleted file mode 100644 index f5c8e5c..0000000 --- a/reachability/load_predicates.m +++ /dev/null @@ -1,118 +0,0 @@ -function pred = load_predicates(plant) -% LOAD_PREDICATES Read predicates.json and resolve rhs_expr into numbers. -% -% Returns: -% pred.constants - struct with T_c0, T_cold0, T_f0, T_standby -% pred.operational_deadbands - struct of predicates (each: A_poly, b_poly, meaning) -% pred.safety_limits - struct of halfspace limits (each: A_poly, b_poly, meaning) -% pred.mode_invariants - struct mapping mode-invariant name to -% its conjoined (A_poly, b_poly, components) -% -% Each halfspace is of the form { x : A_poly * x <= b_poly }. Conjunctions -% (polytopes) stack rows into A_poly / b_poly. -% -% Usage: -% plant = pke_params(); -% pred = load_predicates(plant); -% inv2 = pred.mode_invariants.inv2_holds; -% is_in = all(inv2.A_poly * x <= inv2.b_poly); - - here = fileparts(mfilename('fullpath')); - raw = fileread(fullfile(here, 'predicates.json')); - J = jsondecode(raw); - - % --- Constants used in rhs_expr evaluations --- - T_c0 = plant.T_c0; %#ok - T_f0 = plant.T_f0; %#ok - T_cold0 = plant.T_cold0; %#ok - T_standby_offset_C = J.derived.T_standby_offset_C; - T_standby = T_c0 + T_standby_offset_C; %#ok - - pred.constants = struct( ... - 'T_c0', plant.T_c0, ... - 'T_f0', plant.T_f0, ... - 'T_cold0', plant.T_cold0, ... - 'T_standby', T_standby, ... - 'T_standby_offset_C', T_standby_offset_C, ... - 'T_standby_offset_F', J.derived.T_standby_offset_F, ... - 't_avg_in_range_halfwidth_C', J.derived.t_avg_in_range_halfwidth_C, ... - 'p_above_crit_threshold_n', J.derived.p_above_crit_threshold_n, ... - 'T_fuel_limit_C', J.derived.T_fuel_limit_C, ... - 'T_c_high_trip_C', J.derived.T_c_high_trip_C, ... - 'n_high_trip', J.derived.n_high_trip); - - % --- operational_deadbands --- - pred.operational_deadbands = parse_group(J.operational_deadbands, ... - T_c0, T_f0, T_cold0, T_standby); - - % --- safety_limits --- - pred.safety_limits = parse_group(J.safety_limits, ... - T_c0, T_f0, T_cold0, T_standby); - - % --- mode_invariants: conjunctions of safety_limits entries --- - inv_names = fieldnames(J.mode_invariants); - for k = 1:numel(inv_names) - name = inv_names{k}; - if startsWith(name, '_') || startsWith(name, 'x_'), continue, end - entry = J.mode_invariants.(name); - if ~isstruct(entry) || ~isfield(entry, 'conjunction_of'), continue, end - components = entry.conjunction_of; - if ischar(components), components = {components}; end - A_all = []; - b_all = []; - for i = 1:numel(components) - comp = components{i}; - A_all = [A_all; pred.safety_limits.(comp).A_poly]; %#ok - b_all = [b_all; pred.safety_limits.(comp).b_poly]; %#ok - end - pred.mode_invariants.(name).A_poly = A_all; - pred.mode_invariants.(name).b_poly = b_all; - pred.mode_invariants.(name).meaning = entry.meaning; - pred.mode_invariants.(name).components = components; - end -end - -function group_out = parse_group(group_in, T_c0, T_f0, T_cold0, T_standby) - names = fieldnames(group_in); - group_out = struct(); - for k = 1:numel(names) - name = names{k}; - % MATLAB jsondecode renames "_comment" -> "x_comment" and similar - if startsWith(name, '_') || startsWith(name, 'x_'), continue, end - entry = group_in.(name); - if ~isstruct(entry) || ~isfield(entry, 'halfspaces'), continue, end - hs_list = entry.halfspaces; - if iscell(hs_list) - n_hs = numel(hs_list); - get_hs = @(i) hs_list{i}; - else - n_hs = numel(hs_list); - get_hs = @(i) hs_list(i); - end - A_poly = zeros(n_hs, 10); - b_poly = zeros(n_hs, 1); - for i = 1:n_hs - hs = get_hs(i); - if isfield(hs, 'row') - % Multi-coefficient halfspace: hs.row is a k x 2 matrix - % where each row is [state_index, coeff] (jsondecode from - % JSON array-of-arrays). - coeffs = hs.row; - for r = 1:size(coeffs, 1) - A_poly(i, coeffs(r, 1)) = coeffs(r, 2); - end - else - % Single-coefficient halfspace: hs.state_index + hs.coeff. - A_poly(i, hs.state_index) = hs.coeff; - end - b_poly(i) = evalin_context(hs.rhs_expr, T_c0, T_f0, T_cold0, T_standby); - end - group_out.(name).A_poly = A_poly; - group_out.(name).b_poly = b_poly; - group_out.(name).meaning = entry.meaning; - end -end - -function val = evalin_context(expr, T_c0, T_f0, T_cold0, T_standby) %#ok - val = eval(expr); -end diff --git a/reachability/reach_linear.m b/reachability/reach_linear.m deleted file mode 100644 index 20e84bb..0000000 --- a/reachability/reach_linear.m +++ /dev/null @@ -1,93 +0,0 @@ -function [T, R_lo, R_hi, center] = reach_linear(A_cl, B_w, x0_center, x0_halfwidth, w_lo, w_hi, tspan, dt) -% REACH_LINEAR Interval reach set for dx/dt = A_cl*x + B_w*w, w in [w_lo, w_hi]. -% -% Hand-rolled zonotope propagator specialized to the case where: -% - The initial set is an axis-aligned box around x0_center with -% halfwidths x0_halfwidth. -% - The disturbance w is a scalar interval. -% -% The reach set at time t is: -% X(t) = {expm(A_cl*t) * x0 + integral_0^t expm(A_cl*(t-s))*B_w*w(s) ds : -% x0 in X0, w(s) in W}. -% -% Decompose: -% X(t) = expm(A_cl*t) * X0 (+) S_w(t), -% where S_w(t) is the reachable contribution of the disturbance. -% Both parts are zonotopes; we over-approximate each by a box at every -% time step via the matrix/integral absolute sum (interval hull). -% -% Inputs: -% A_cl - n x n closed-loop state matrix (Hurwitz assumed) -% B_w - n x 1 disturbance vector -% x0_center - n x 1 center of initial set -% x0_halfwidth - n x 1 halfwidths of initial set (>=0) -% w_lo, w_hi - scalar disturbance bounds -% tspan - [t0, tf] -% dt - time step for exporting the reach tube (s) -% -% Outputs: -% T - M x 1 time grid -% R_lo - n x M lower bounds of reach set -% R_hi - n x M upper bounds of reach set -% center- n x M center trajectory (nominal, w = w_mid) - - n = size(A_cl, 1); - assert(size(A_cl,2) == n, 'A_cl must be square'); - assert(numel(x0_center) == n && numel(x0_halfwidth) == n, 'x0 dim mismatch'); - x0_center = x0_center(:); - x0_halfwidth = x0_halfwidth(:); - B_w = B_w(:); - - w_mid = 0.5 * (w_lo + w_hi); - w_half = 0.5 * (w_hi - w_lo); - - T = (tspan(1):dt:tspan(2)).'; - M = numel(T); - - R_lo = zeros(n, M); - R_hi = zeros(n, M); - center = zeros(n, M); - - % Disturbance contribution: - % S_center(t) = int_0^t expm(A*(t-s))*B_w*w_mid ds (signed, centerline) - % S_half(t) = int_0^t |expm(A*(t-s))*B_w| * w_half ds (elementwise, halfwidth) - % - % The halfwidth uses elementwise |.| because under an axis-aligned - % box over-approximation, |linear-map| contributes absolute values. - % This is the interval-arithmetic analog of propagating a zonotope. - S_center = zeros(n, 1); - S_half = zeros(n, 1); - - x_center_now = x0_center; - A_step = expm(A_cl * dt); % state transition per step - A_abs_step = abs(A_step); % for interval-halfwidth - - % Integrated input gain per step: - % G_step = int_0^dt expm(A*s)*B_w ds - % Augmented-matrix trick: upper-right block of expm([A B_w; 0 0]*dt). - M_aug = expm([A_cl, B_w; zeros(1, n+1)] * dt); - G_step = M_aug(1:n, n+1); % n x 1 - G_abs_step = abs(G_step); - - halfwidth_now = x0_halfwidth; - - % Initial state - center(:, 1) = x_center_now + S_center; - R_lo(:, 1) = center(:, 1) - halfwidth_now - S_half; - R_hi(:, 1) = center(:, 1) + halfwidth_now + S_half; - - for k = 2:M - % State piece - x_center_now = A_step * x_center_now; - halfwidth_now = A_abs_step * halfwidth_now; - - % Disturbance piece - S_center = A_step * S_center + G_step * w_mid; - S_half = A_abs_step * S_half + G_abs_step * w_half; - - center(:, k) = x_center_now + S_center; - R_lo(:, k) = center(:, k) - halfwidth_now - S_half; - R_hi(:, k) = center(:, k) + halfwidth_now + S_half; - end - -end diff --git a/reachability/reach_operation.m b/reachability/reach_operation.m deleted file mode 100644 index 0dd854f..0000000 --- a/reachability/reach_operation.m +++ /dev/null @@ -1,203 +0,0 @@ -%% reach_operation.m — linear reach set for operation mode (LQR closed-loop) -% -% *** SOUNDNESS STATUS: APPROXIMATE, NOT SOUND. *** -% -% This file computes a reach tube for the *linearized* closed-loop -% system (A_cl = A - BK around x_op) under bounded Q_sg. The tube -% itself is a sound over-approximation of the LINEAR model's reach -% set — it uses conservative box hulls and elementwise-absolute-value -% matrix propagation. But the LINEAR MODEL is only an approximation -% of the real nonlinear plant (pke_th_rhs.m), so the result is not a -% sound reach tube for the actual plant. -% -% To upgrade to a sound result, pick one: -% (a) Nonlinear reach directly (CORA's nonlinearSys, JuliaReach's -% BlackBoxContinuousSystem). Expensive, honest. -% (b) Linear reach + Taylor-remainder inflation: compute an upper -% bound on ||f_nonlinear(x,u) - (A x + B u)|| over the reach -% set and inflate the linear tube by that bound. Cheaper, -% requires a Hessian-norm estimate. -% Tracked as a thesis-blocking todo. For now, the 5-orders-of-margin -% buffer (|dT_c| ~ 0.03 K vs safety band 5 K) gives us a lot of room -% to absorb linearization error, but that's not a proof. -% -% Compute a reach-tube over-approximation starting from a box around -% x_op, under LQR feedback, with Q_sg in a specified interval. Check -% that T_avg stays inside the t_avg_in_range predicate for all t in -% [0, T_final]. -% -% This is the *continuous-mode obligation* for q_operation: -% X_entry := { x : |x - x_op| <= delta_entry } -% W := [Q_min, Q_max] -% X_safe := { x : |T_c - T_c0| <= delta_safe } -% Obligation: ReachTube(X_entry, W, [0, T_final]) subset X_safe. -% -% If this passes, we've discharged the Thrust-3 verification for one -% continuous mode at the level the thesis calls for. - -clear; clc; close all; -addpath('../plant-model', '../plant-model/controllers'); - -plant = pke_params(); -x_op = pke_initial_conditions(plant); -pred = load_predicates(plant); % single source of truth for predicate bands - -%% ===== Closed-loop linearization ===== -[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); - -% LQR gain using the same Q, R as ctrl_operation_lqr.m. Keep these in -% sync if you retune there. -Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1]); -R_lqr = 1e6; -try - K = lqr(A, B, Q_lqr, R_lqr); -catch - [~, ~, K] = icare(A, B, Q_lqr, R_lqr); -end - -A_cl = A - B*K; - -fprintf('\n=== Closed-loop spectrum (A - BK) ===\n'); -eigs_cl = eig(A_cl); -fprintf(' max Re(eig) = %.3e\n', max(real(eigs_cl))); -fprintf(' min Re(eig) = %.3e\n', min(real(eigs_cl))); -assert(all(real(eigs_cl) < -1e-8), 'A_cl not Hurwitz — gain tuning issue'); - -%% ===== Define sets ===== -% X_entry: 1% box on n, 0.1% boxes on precursors (their magnitudes are -% huge due to 1/Lambda), 0.1 K on fuel, 0.1 K on coolant, 0.1 K on cold leg. -% This represents "we entered operation mode near the steady state". -delta_entry = [0.01 * x_op(1); % n - 0.001 * abs(x_op(2:7)); % C1..C6 - 0.1; % T_f [C] - 0.1; % T_c [C] - 0.1]; % T_cold [C] - -% Disturbance: Q_sg = [0.85*P0, 1.00*P0] -> this captures up to a 15% -% down-step load demand (realistic load-follow envelope). -Q_nom = plant.P0; -Q_min = 0.85 * plant.P0; -Q_max = 1.00 * plant.P0; -dQ_lo = Q_min - Q_nom; % -0.15 * P0 -dQ_hi = Q_max - Q_nom; % 0 - -% X_safe is inv2_holds — the operation-mode safety envelope. Each row of -% inv2.A_poly is a hard safety limit (fuel centerline, T_c high/low trip, -% n high/low, cold-leg subcooling). We check reach-tube containment -% halfspace-by-halfspace so failure modes are attributable. -inv2 = pred.mode_invariants.inv2_holds; -% Keep the old +-halfwidth report too, for continuity with the -% operational-deadband framing. -delta_safe_Tc = pred.constants.t_avg_in_range_halfwidth_C; % [C] - -%% ===== Reach set ===== -% Propagate in deviation coordinates: dx = x - x_op. -tspan = [0, 600]; -dt = 0.5; - -[T, R_lo, R_hi, C] = reach_linear(A_cl, B_w, zeros(10,1), delta_entry, dQ_lo, dQ_hi, tspan, dt); - -% Translate back to absolute coordinates for reporting -Xabs_lo = R_lo + x_op; -Xabs_hi = R_hi + x_op; -Cabs = C + x_op; - -%% ===== Safety check ===== -T_c_lo = Xabs_lo(9, :); -T_c_hi = Xabs_hi(9, :); - -violation_mask = (T_c_hi > plant.T_c0 + delta_safe_Tc) | ... - (T_c_lo < plant.T_c0 - delta_safe_Tc); -fprintf('\n=== Operation-mode reach-set safety ===\n'); -fprintf(' Horizon = [%g, %g] s\n', tspan(1), tspan(2)); -fprintf(' Entry box T_c [C] = [%.3f, %.3f] (x_op +/- %.1f C)\n', ... - x_op(9) - delta_entry(9), x_op(9) + delta_entry(9), delta_entry(9)); -fprintf(' Disturbance Q_sg = [%.3f, %.3f] MW\n', Q_min/1e6, Q_max/1e6); -fprintf(' Safe band on T_c = x_op(T_c0) +/- %.1f C -> [%.3f, %.3f]\n', ... - delta_safe_Tc, plant.T_c0 - delta_safe_Tc, plant.T_c0 + delta_safe_Tc); -fprintf(' Reach T_c envelope = [%.3f, %.3f]\n', min(T_c_lo), max(T_c_hi)); -if any(violation_mask) - t_first = T(find(violation_mask, 1)); - fprintf(' *** SAFETY VIOLATED at t = %.2f s ***\n', t_first); -else - fprintf(' OK: reach set stays inside the safe band.\n'); -end - -%% Hard safety-limit check (inv2_holds halfspace-by-halfspace) -% For each row a_k of inv2.A_poly with threshold b_k, check whether -% max over reach tube of a_k * x stays <= b_k. The reach tube upper -% envelope is Xabs_hi; lower envelope is Xabs_lo. We evaluate -% max(a_k * x) using Xabs_hi where a_k > 0, Xabs_lo where a_k < 0. -fprintf('\n=== Operation-mode reach vs inv2_holds safety limits ===\n'); -A_inv = inv2.A_poly; b_inv = inv2.b_poly; -comps = inv2.components; -for k = 1:size(A_inv, 1) - a = A_inv(k, :).'; - % envelope maximum of a' * x across the reach tube - x_envelope = Xabs_hi .* (a > 0) + Xabs_lo .* (a < 0); % 10 x M - max_ax = max(a.' * x_envelope); - margin = b_inv(k) - max_ax; - status = 'OK'; - if margin < 0, status = '*** VIOLATED ***'; end - fprintf(' [%s] a''x <= %.3f | max a''x = %.3f | margin = %+.3f %s\n', ... - comps{k}, b_inv(k), max_ax, margin, status); -end - -%% Per-state reach-set growth diagnostic (final time vs initial) -state_names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'}; -fprintf('\n=== Reach-set width at t=0 vs t=T_final ===\n'); -fprintf(' %-7s %-14s %-14s %-8s\n', 'state', 'init halfwidth', 'final halfwidth', 'ratio'); -for i = 1:10 - hi = 0.5 * (R_hi(i, 1) - R_lo(i, 1)); - hf = 0.5 * (R_hi(i, end) - R_lo(i, end)); - fprintf(' %-7s %-14.4e %-14.4e %-8.2f\n', state_names{i}, hi, hf, hf/max(hi,eps)); -end - -%% ===== Plots ===== -figdir = fullfile('..', 'docs', 'figures'); -if ~exist(figdir, 'dir'), mkdir(figdir); end -CtoF = @(T) T*9/5 + 32; - -% Two-panel plot: wide view with safety band, zoom view showing actual tube. -figure('Position', [100 80 1400 500], 'Name', 'Reach tube: T_c'); - -subplot(1,2,1); -fill([T; flipud(T)], CtoF([T_c_hi.'; flipud(T_c_lo.')]), [1.0 0.85 0.85], ... - 'EdgeColor', 'none'); hold on; -plot(T, CtoF(Cabs(9, :)), 'r-', 'LineWidth', 1.2); -yline(CtoF(plant.T_c0 + delta_safe_Tc), 'k--', 'LineWidth', 1.0); -yline(CtoF(plant.T_c0 - delta_safe_Tc), 'k--', 'LineWidth', 1.0); -yline(CtoF(plant.T_c0), 'k:', 'LineWidth', 1.0); -xlabel('Time [s]'); ylabel('T_{avg} [F]'); grid on; -title('Safety-band view'); -legend('reach tube', 'nominal', 'safety +/- 5 C', 'Location', 'best'); - -subplot(1,2,2); -Tc_dev_lo = T_c_lo.' - plant.T_c0; % M x 1, deviation in K -Tc_dev_hi = T_c_hi.' - plant.T_c0; -fill([T; flipud(T)], [Tc_dev_hi; flipud(Tc_dev_lo)], [1.0 0.85 0.85], ... - 'EdgeColor', 'none'); hold on; -plot(T, Cabs(9, :).' - plant.T_c0, 'r-', 'LineWidth', 1.2); -yline(0, 'k:', 'LineWidth', 1.0); -xlabel('Time [s]'); ylabel('T_{avg} - T_{c0} [K]'); grid on; -max_dev = max(abs([Tc_dev_lo; Tc_dev_hi])); -title(sprintf('Zoomed: max |dT_c| = %.3e K', max_dev)); - -sgtitle(sprintf('Operation-mode reach tube, LQR, Q_{sg} in [%.0f%%, %.0f%%] P_0', ... - 100*Q_min/Q_nom, 100*Q_max/Q_nom)); -exportgraphics(gcf, fullfile(figdir, 'reach_operation_Tc.png'), 'Resolution', 150); - -figure('Position', [100 80 1100 500], 'Name', 'Reach tube: n'); -fill([T; flipud(T)], [R_hi(1,:).' + x_op(1); flipud(R_lo(1,:).' + x_op(1))], ... - [0.85 0.85 1.0], 'EdgeColor', 'none'); hold on; -plot(T, Cabs(1, :), 'b-', 'LineWidth', 1.2); -xlabel('Time [s]'); ylabel('n'); grid on; -title('Operation mode reach tube on normalized power'); -legend('reach tube', 'nominal', 'Location', 'best'); -exportgraphics(gcf, fullfile(figdir, 'reach_operation_n.png'), 'Resolution', 150); - -save(fullfile('.', 'reach_operation_result.mat'), ... - 'T', 'R_lo', 'R_hi', 'C', 'Xabs_lo', 'Xabs_hi', 'Cabs', ... - 'K', 'A_cl', 'x_op', 'delta_entry', 'Q_min', 'Q_max', 'delta_safe_Tc', '-v7'); - -fprintf('\nSaved reach result to ./reach_operation_result.mat\n');