PWR-HYBRID-3/code/scripts/sim/sim_sanity.jl
Dane Sabo 2bbb1871cc refactor: scripts subdivision + TOML configs + results/ split + presentation outline
Architecture restructure from morning review:

1. code/scripts/ subdivided into sim/, reach/, barrier/, plot/.
   Easier nav; `barrier/` is the natural place for SOS scale-up scripts.
2. Heatup PJ reach variants consolidated behind TOML configs.
   reach_heatup_pj.jl now takes `--config path/to/config.toml`;
   configs/heatup/baseline.toml (wide entry, from predicates.json) and
   configs/heatup/tight.toml (narrow entry, reproduces all-6-halfspaces
   discharged result). Old reach_heatup_pj_tight.jl and
   reach_heatup_pj_tight_full.jl deleted (superseded).
3. Reach output .mat files moved from reachability/ to results/.
   reachability/ now = specs + docs; results/ = ephemeral outputs
   (gitignored *.mat). README added.
4. OVERNIGHT_NOTES.md archived to claude_memory/2026-04-20-21-overnight-
   session-summary.md (date range in the filename makes the history clearer).

All include() / Pkg.activate() paths in scripts updated for the new
depth. Smoke tests pass (reach_operation.jl generates its .mat in
the new results/ location; sim_sanity.jl matches MATLAB).

Presentation outline for the 20-min prelim talk landed in
presentations/prelim-presentation/outline.md. 14-slide assertion-
evidence format targeting OT-informed cybersecurity audience. Each
slide: one declarative assertion + one figure. Outline includes
which figures already exist and which need to be created, timing
checkpoints, cybersecurity angle to emphasize, and Q&A prep.

New config configs/heatup/with_steam_dump.toml + its companion
scripts/reach/reach_heatup_pj_sd.jl (12-state RHS with Q_sg as an
augmented bounded parameter x[10] and time as x[11]). Kicks off
point 3 from morning review.

Next up: scram X_entry expansion (morning point 2) — LOCA scenario
+ union of mode reach envelopes.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
2026-04-21 20:24:48 -04:00

41 lines
1.4 KiB
Julia

#!/usr/bin/env julia
#
# sim_sanity.jl — verify the Julia port matches the MATLAB result.
#
# Reproduces main.m Run 2 (ctrl_operation under 100% -> 80% Q_sg step)
# and prints the final state, which should agree with MATLAB to ~1e-3.
using Pkg
Pkg.activate(joinpath(@__DIR__, "..", ".."))
using OrdinaryDiffEq
include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl"))
include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl"))
include(joinpath(@__DIR__, "..", "..", "controllers", "controllers.jl"))
plant = pke_params()
x0 = pke_initial_conditions(plant)
Q_sg = t -> plant.P0 * (1.0 - 0.2 * (t >= 30))
ref = (; T_avg = plant.T_c0)
function rhs!(dx, x, p, t)
u = ctrl_operation(t, x, plant, ref)
pke_th_rhs!(dx, x, t, plant, Q_sg, u)
end
prob = ODEProblem(rhs!, x0, (0.0, 600.0))
sol = solve(prob, Rodas5(); reltol=1e-8, abstol=1e-10)
xf = sol.u[end]
CtoF(T) = T * 9/5 + 32
println("\n=== Julia port sanity — ctrl_operation under 100% -> 80% Q_sg step ===")
println(" Final t = ", sol.t[end])
println(" n = $(round(xf[1]; digits=4)) (expect ~0.800)")
println(" T_f = $(round(CtoF(xf[8]); digits=2)) F (expect ~616.6)")
println(" T_avg = $(round(CtoF(xf[9]); digits=2)) F (expect ~587.8)")
println(" T_cold = $(round(CtoF(xf[10]); digits=2)) F (expect ~561.4)")
u_final = ctrl_operation(sol.t[end], xf, plant, ref)
println(" u = $(round(u_final/plant.beta; digits=4)) \$ (expect ~-0.0068)")