PWR-HYBRID-3/code/CLAUDE.md

# 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_states.jl            canonical ICs for every DRC mode (single source)
    pke_th_rhs.jl            dynamics f(t, x, plant, Q_sg, u)
    pke_linearize.jl         numerical (A, B, B_w) Jacobians
    pke_solver.jl            closed-loop OrdinaryDiffEq driver
    plot_pke_results.jl      4-panel results plot (Plots.jl)
    reach_linear.jl          hand-rolled box reach propagator
    load_predicates.jl       reads reachability/predicates.json
    sos_barrier.jl           solve_sos_barrier_2d helper (SOS feasibility +
                             ε-slack pattern; used by barrier_sos_2d{,_shutdown}.jl)
  controllers/
    controllers.jl           ctrl_null, shutdown, heatup, operation,
                             operation_lqr factory, scram
    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.
- **Single-point ICs live in `pke_states(plant; predicates=...)`.** Returns
  a NamedTuple with `T_standby`, `op`, `shutdown`, `heatup` — every script
  that needs an IC for forward sim or equilibrium-finding pulls from here.
  X_entry *polytopes* (for reach) live in `predicates.json` and are loaded
  separately by reach scripts; don't conflate the two.
- **SOS barrier programs need ε-slack + scale cap.** The naive feasibility
  formulation silently returns `B ≡ 0`. `solve_sos_barrier_2d` in
  `src/sos_barrier.jl` bakes in the (ε > 0, capped) pattern. Use it for
  any new SOS barrier work.
- **Plant parameters live in `pke_params()` as a NamedTuple.** For
  reach analysis (TMJets) we duplicate them as `const` globals in
  each reach script — `@taylorize` needs compile-time constants, not
  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.