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:

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

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

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.