danesabo/PWR-HYBRID-3
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

prompt-jump model + app v2 + overnight journal entry (in progress)

Browse Source 
Singular-perturbation reduction of the PKE+T/H system: set dn/dt=0,
solve algebraically n = Λ·Σλ_i·C_i / (β-ρ). State drops 10 -> 9 (no
n), removes Λ⁻¹ stiffness. Validated against full state on the heatup
scenario:

  t [s]    |Δn|/n_full   T_c err [K]
  60       3.7e-5        4e-6
  300      3.8e-4        1.9e-4
  1200     1.0e-3        2.2e-3
  3000     5.0e-4        7.2e-3

Maximum relative error 0.1% on n, peak 7 mK on temperatures over
50 minutes.  PJ approximation is excellent for slow heatup transients
(sub-prompt-critical regime).

Files:
  - code/src/pke_th_rhs_pj.jl: reduced 9-state RHS
  - code/scripts/validate_pj.jl: side-by-side sim
  - code/scripts/reach_heatup_pj.jl: TMJets reach with PJ model
    (probing T = 60, 300, 1800, 5400 s)

App v2 (Pluto):
  - §9b: live ingestion of reach_operation_result.mat with per-
    halfspace margins computed from JSON-defined inv2_holds.
  - §9c: 2D projection chooser (n, T_f, T_c, T_cold) with reach
    tube envelope overlay.
  - §9d: PJ heatup reach summary (placeholder until first run lands).

Journal:
  - Added 2026-04-20-overnight-prompt-jump.tex with PJ derivation,
    validation table, soundness ledger update.  apass markers for
    the in-progress reach results.

This commit captures state mid-run; next commit will add the
populated reach results once TMJets returns.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This commit is contained in:
Dane Sabo 2026-04-20 22:45:24 -04:00
parent 5acaa5553d
commit 645f2d8d27
8 changed files with 825 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
app/Project.toml
View File

@ -2,6 +2,7 @@ authors = ["Dane Sabo <yourstruly@danesabo.com>"]
[deps]
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

169
app/predicate_explorer.jl
View File

@ -26,6 +26,7 @@ begin
using JSON
using PlutoUI
using Plots
using MAT
gr()
end
@ -312,6 +313,165 @@ md"""
| `q_shutdown` | (TBD) | ❌ not started (trivial — constant u) | — |
"""
# ╔═╡ a74106fb-d501-428a-a02e-41f55b49b3da
md"""
## 9b · Live reach result ingestion
If `code/scripts/reach_operation.jl` has been run, `reach_operation_result.mat`
exists in `../reachability/`. Load it and show per-halfspace margins live —
no more hand-maintained traceability table.
"""
# ╔═╡ 37d6b212-9f00-4684-9f91-50c7e17cbd62
begin
reach_mat_path = joinpath(@__DIR__, "..", "reachability", "reach_operation_result.mat")
have_reach_mat = isfile(reach_mat_path)
if have_reach_mat
reach_mat = matread(reach_mat_path)
md"""**Loaded** `reach_operation_result.mat` — keys: $(join(sort(collect(keys(reach_mat))), ", "))."""
else
md"""_(no `reach_operation_result.mat` found — run `cd code && julia --project=. scripts/reach_operation.jl`)_"""
end
end
# ╔═╡ df53126a-1efa-4cc7-87f4-be4ff0808513
begin
inv2_haspath = haskey(pred_raw["mode_invariants"], "inv2_holds")
if have_reach_mat && inv2_haspath
# Reconstruct per-halfspace margins from the loaded reach tube.
T_op_const = 308.349
T_f0_const = 328.349
Tcold0_const = 290.0
Tstandby_local = T_op_const + pred_raw["derived"]["T_standby_offset_C"]
eval_local(expr) = let T_c0=T_op_const, T_f0=T_f0_const, T_cold0=Tcold0_const, T_standby=Tstandby_local
eval(Meta.parse(expr))
end
# Build A_inv, b_inv from JSON for inv2_holds (same logic as load_predicates).
inv2_entry = pred_raw["mode_invariants"]["inv2_holds"]
components = inv2_entry["conjunction_of"]
A_inv = zeros(0, 10)
b_inv = zeros(0)
labels = String[]
for comp in components
entry = pred_raw["safety_limits"][comp]
for hs in entry["halfspaces"]
row = zeros(1, 10)
if haskey(hs, "row")
for pair in hs["row"]
row[1, Int(pair[1])] = Float64(pair[2])
end
else
row[1, Int(hs["state_index"])] = Float64(hs["coeff"])
end
A_inv = vcat(A_inv, row)
b = eval_local(hs["rhs_expr"])
b_inv = vcat(b_inv, b)
push!(labels, comp)
end
end
X_lo_mat = reach_mat["X_lo"]
X_hi_mat = reach_mat["X_hi"]
rows = String[]
push!(rows, "| Safety halfspace | Limit | Reach max | Margin | |")
push!(rows, "|---|---:|---:|---:|---|")
for k in 1:size(A_inv, 1)
a = A_inv[k, :]
env = (a .> 0) .* X_hi_mat .+ (a .< 0) .* X_lo_mat
zero_mask = (a .== 0)
env[zero_mask, :] .= X_hi_mat[zero_mask, :]
max_ax = maximum(a' * env)
margin = b_inv[k] - max_ax
badge = margin >= 0 ? "✅" : "❌"
push!(rows, "| `$(labels[k])` | $(round(b_inv[k]; digits=3)) | $(round(max_ax; digits=3)) | $(round(margin; digits=3)) | $badge |")
end
Markdown.parse(join(rows, "\n"))
else
md"_(reach result or inv2_holds unavailable; cannot render live margins)_"
end
end
# ╔═╡ 9bd8d609-ffa4-46a7-8f0f-d093d87e45e9
md"""
## 9c · 2D projection chooser
Pick any two state coordinates; render the operating polytope inside the
safety limits, with the reach tube envelope overlaid (if loaded).
"""
# ╔═╡ e4d56600-7a46-46a1-a7c9-f65a5db95f66
@bind proj_x Select(["n", "T_f", "T_c", "T_cold"], default="T_c")
# ╔═╡ 4fdf70de-3f58-4a38-917b-62b3689ce534
@bind proj_y Select(["n", "T_f", "T_c", "T_cold"], default="n")
# ╔═╡ 37251c95-dd46-402f-9579-d4ede2c1e5ff
begin
state_idx = Dict("n"=>1, "T_f"=>8, "T_c"=>9, "T_cold"=>10)
state_units = Dict("n"=>"", "T_f"=>" [°C]", "T_c"=>" [°C]", "T_cold"=>" [°C]")
px_idx = state_idx[proj_x]
py_idx = state_idx[proj_y]
# Determine reasonable axis ranges from operating point + 5% margins.
x_op_demo = [1.0,
173.4, 467.0, 114.8, 85.3, 6.56, 0.91,
328.35, 308.35, 290.0]
px_op = x_op_demo[px_idx]
py_op = x_op_demo[py_idx]
px_span = max(0.2 * abs(px_op), 5.0)
py_span = max(0.2 * abs(py_op), 5.0)
p_proj = plot(xlim=(px_op - px_span, px_op + px_span),
ylim=(py_op - py_span, py_op + py_span),
xlabel=proj_x * state_units[proj_x],
ylabel=proj_y * state_units[proj_y],
title="Projection: $proj_y vs $proj_x",
size=(700, 450), legend=:topleft)
# Overlay reach tube envelope as a rectangle.
if have_reach_mat
X_lo_local = reach_mat["X_lo"]
X_hi_local = reach_mat["X_hi"]
rx_lo = minimum(X_lo_local[px_idx, :])
rx_hi = maximum(X_hi_local[px_idx, :])
ry_lo = minimum(X_lo_local[py_idx, :])
ry_hi = maximum(X_hi_local[py_idx, :])
plot!(p_proj, [rx_lo, rx_hi, rx_hi, rx_lo, rx_lo],
[ry_lo, ry_lo, ry_hi, ry_hi, ry_lo],
lw=2, color=:red, label="reach tube envelope")
end
# Operating point.
scatter!(p_proj, [px_op], [py_op], color=:black, markersize=8,
label="x_op")
p_proj
end
# ╔═╡ 4d8459c1-a937-4b53-9eca-975261d6a2cd
md"""
## 9d · Prompt-jump heatup reach (if available)
If `code/scripts/reach_heatup_pj.jl` has been run, the latest result
shows up here as a T_c reach envelope vs time.
This is the singular-perturbation reduction that lets nonlinear reach
extend past the ~10 s prompt-neutron stiffness wall.
"""
# ╔═╡ 45e8962a-873e-48d3-aac4-70ea0cf36c97
begin
pj_mat_path = joinpath(@__DIR__, "..", "reachability", "reach_heatup_pj_result.mat")
if isfile(pj_mat_path)
pj_mat = matread(pj_mat_path)
md"""**Loaded** `reach_heatup_pj_result.mat` — pending plot rendering."""
else
md"""_(no `reach_heatup_pj_result.mat` yet — run `cd code && julia --project=. scripts/reach_heatup_pj.jl`)_"""
end
end
# ╔═╡ de922a70-b653-4b2c-bb13-ccc14b14ac91
md"""
## 10 · Notes for the next pass
@ -357,4 +517,13 @@ md"""
# ╟─b7fdf92d-948b-41f2-8516-dd8fc20c2efe
# ╟─0560c66a-c12f-4f15-bc49-2c9c8164abd1
# ╟─1f6abf31-3618-46a5-9f57-4cb3c8022f89
# ╟─a74106fb-d501-428a-a02e-41f55b49b3da
# ╠═37d6b212-9f00-4684-9f91-50c7e17cbd62
# ╠═df53126a-1efa-4cc7-87f4-be4ff0808513
# ╟─9bd8d609-ffa4-46a7-8f0f-d093d87e45e9
# ╟─e4d56600-7a46-46a1-a7c9-f65a5db95f66
# ╟─4fdf70de-3f58-4a38-917b-62b3689ce534
# ╠═37251c95-dd46-402f-9579-d4ede2c1e5ff
# ╟─4d8459c1-a937-4b53-9eca-975261d6a2cd
# ╠═45e8962a-873e-48d3-aac4-70ea0cf36c97
# ╟─de922a70-b653-4b2c-bb13-ccc14b14ac91

218
code/scripts/reach_heatup_pj.jl Normal file
View File

@ -0,0 +1,218 @@
#!/usr/bin/env julia
#
# reach_heatup_pj.jl — nonlinear reach on heatup, prompt-jump model.
#
# Reduced from 10-state to 9-state (n is algebraic). Removes the
# Λ⁻¹ stiffness that capped the full-state reach at ~10 s. We push
# horizons up: 60 s, 300 s, 1800 s, 5400 s, full T_max = 18000 s (5 hr).
#
# State (10D with augmented time):
# x[1..6] = C_1..C_6 (delayed-neutron precursors)
# x[7] = T_f
# x[8] = T_c
# x[9] = T_cold
# x[10] = t (augmented time, dt/dt = 1)
#
# n is algebraic: n = Λ · Σ λ_i C_i / (β - ρ), ρ = K_p · (T_ref - T_c).
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using LinearAlgebra
using ReachabilityAnalysis, LazySets
using JSON
using MAT
# --- Inlined plant constants (must match pke_params) ---
const LAMBDA = 1e-4
const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 =
0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273
const BETA = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 =
0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01
const P0 = 1e9
const M_F, C_F, M_C, C_C, HA, W_M, M_SG =
50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0
# Note: feedback-linearization in ctrl_heatup_unsat cancels the alpha
# terms exactly, so under closed-loop the effective rho is just Kp·e.
# We don't need ALPHA_F, ALPHA_C in the reach RHS as a result.
const T_COLD0 = 290.0
const DT_CORE = P0 / (W_M * C_C)
const T_HOT0 = T_COLD0 + DT_CORE
const T_C0 = (T_HOT0 + T_COLD0) / 2
const T_F0 = T_C0 + P0 / HA
const T_STANDBY = T_C0 - 33.333333
const RAMP_RATE_CS = 28.0 / 3600
const KP_HEATUP = 1e-4
# --- Taylorized PJ heatup RHS, 10D with augmented time ---
@taylorize function rhs_heatup_pj_taylor!(dx, x, p, t)
# rho_total under closed-loop feedback linearization = Kp · (T_ref - T_c).
rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8])
# Algebraic prompt-jump n.
sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] +
LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6]
denom = BETA - rho
n = LAMBDA * sum_lam_C / denom
inv_factor = sum_lam_C / denom
# Precursor balance under PJ.
dx[1] = BETA_1 * inv_factor - LAM_1 * x[1]
dx[2] = BETA_2 * inv_factor - LAM_2 * x[2]
dx[3] = BETA_3 * inv_factor - LAM_3 * x[3]
dx[4] = BETA_4 * inv_factor - LAM_4 * x[4]
dx[5] = BETA_5 * inv_factor - LAM_5 * x[5]
dx[6] = BETA_6 * inv_factor - LAM_6 * x[6]
# Thermal — n appears algebraically in fuel eq.
dx[7] = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F)
dx[8] = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C)
dx[9] = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C)
# Augmented time.
dx[10] = one(x[1])
return nothing
end
# --- Build X_entry (PJ form: no n) from predicates.json ---
pred_path = joinpath(@__DIR__, "..", "..", "reachability", "predicates.json")
pred_raw = JSON.parsefile(pred_path)
entry = pred_raw["mode_boundaries"]["q_heatup"]["X_entry_polytope"]
n_lo, n_hi = entry["n_range"]
T_f_lo, T_f_hi = entry["T_f_range_C"]
T_c_lo, T_c_hi = entry["T_c_range_C"]
T_cold_lo, T_cold_hi = entry["T_cold_range_C"]
n_mid = 0.5 * (n_lo + n_hi)
C_mid_1 = (BETA_1 / (LAM_1 * LAMBDA)) * n_mid
C_mid_2 = (BETA_2 / (LAM_2 * LAMBDA)) * n_mid
C_mid_3 = (BETA_3 / (LAM_3 * LAMBDA)) * n_mid
C_mid_4 = (BETA_4 / (LAM_4 * LAMBDA)) * n_mid
C_mid_5 = (BETA_5 / (LAM_5 * LAMBDA)) * n_mid
C_mid_6 = (BETA_6 / (LAM_6 * LAMBDA)) * n_mid
# C_i scale linearly with n; sweep across the n_lo..n_hi band.
x_lo = [C_mid_1 * (n_lo / n_mid); C_mid_2 * (n_lo / n_mid);
C_mid_3 * (n_lo / n_mid); C_mid_4 * (n_lo / n_mid);
C_mid_5 * (n_lo / n_mid); C_mid_6 * (n_lo / n_mid);
T_f_lo; T_c_lo; T_cold_lo;
0.0]
x_hi = [C_mid_1 * (n_hi / n_mid); C_mid_2 * (n_hi / n_mid);
C_mid_3 * (n_hi / n_mid); C_mid_4 * (n_hi / n_mid);
C_mid_5 * (n_hi / n_mid); C_mid_6 * (n_hi / n_mid);
T_f_hi; T_c_hi; T_cold_hi;
0.0]
X0 = Hyperrectangle(low=x_lo, high=x_hi)
println("\n=== Nonlinear heatup reach, prompt-jump model ===")
println(" X_entry (n-implied range): n ∈ [$(n_lo), $(n_hi)]")
println(" T_c ∈ [$(T_c_lo), $(T_c_hi)] °C")
T_RAMP_END = (T_C0 - T_STANDBY) / RAMP_RATE_CS
println(" T_ramp_end = $(round(T_RAMP_END; digits=0)) s ($(round(T_RAMP_END/60; digits=1)) min)")
println(" Probing horizons up to T_max(heatup) = 18000 s (5 hr).")
# Probe at increasing horizons. Stop early if any probe fails.
probe_horizons = (60.0, 300.0, 1800.0, 5400.0)
results = Dict{Float64, Any}()
for T_probe in probe_horizons
println("\n--- Probe T = $T_probe s ($(round(T_probe/60; digits=1)) min) ---")
sys = BlackBoxContinuousSystem(rhs_heatup_pj_taylor!, 10)
prob = InitialValueProblem(sys, X0)
try
alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000)
t_start = time()
sol = solve(prob; T=T_probe, alg=alg)
elapsed = time() - t_start
flow = flowpipe(sol)
n_sets = length(flow)
println(" TMJets: $n_sets reach-sets, wall-time $(round(elapsed; digits=1)) s")
flow_hr = overapproximate(flow, Hyperrectangle)
# Envelope at the FINAL set.
Tc_lo_env = minimum(low(set(R), 8) for R in flow_hr)
Tc_hi_env = maximum(high(set(R), 8) for R in flow_hr)
Tf_lo_env = minimum(low(set(R), 7) for R in flow_hr)
Tf_hi_env = maximum(high(set(R), 7) for R in flow_hr)
Tcold_lo = minimum(low(set(R), 9) for R in flow_hr)
Tcold_hi = maximum(high(set(R), 9) for R in flow_hr)
# Reconstruct n envelope at each step from C and T_c via PJ formula.
n_env_lo = Inf
n_env_hi = -Inf
for R in flow_hr
s = set(R)
sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) +
LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6)
sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) +
LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6)
# rho range: ρ = Kp*(T_ref - T_c). T_ref bounded by [T_STANDBY, T_C0],
# T_c bounded by current envelope.
rho_lo = KP_HEATUP * (T_STANDBY - high(s, 8))
rho_hi = KP_HEATUP * (T_C0 - low(s, 8))
denom_lo = BETA - rho_hi # smaller denom => larger n
denom_hi = BETA - rho_lo
if denom_lo > 0
n_lo_local = LAMBDA * sumLC_lo / denom_hi
n_hi_local = LAMBDA * sumLC_hi / denom_lo
n_env_lo = min(n_env_lo, n_lo_local)
n_env_hi = max(n_env_hi, n_hi_local)
end
end
println(" n envelope (reconstructed): [$(round(n_env_lo; sigdigits=4)), $(round(n_env_hi; sigdigits=4))]")
println(" T_f envelope: [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C")
println(" T_c envelope: [$(round(Tc_lo_env; digits=2)), $(round(Tc_hi_env; digits=2))] °C")
println(" T_cold envelope: [$(round(Tcold_lo; digits=2)), $(round(Tcold_hi; digits=2))] °C")
results[T_probe] = (status="OK", n_sets=n_sets, elapsed=elapsed,
Tc=(Tc_lo_env, Tc_hi_env),
Tf=(Tf_lo_env, Tf_hi_env),
Tcold=(Tcold_lo, Tcold_hi),
n=(n_env_lo, n_env_hi))
catch err
msg = sprint(showerror, err)
println(" FAILED: ", first(msg, 300))
results[T_probe] = (status="FAILED", err=first(msg, 300))
break
end
end
println("\n=== Summary ===")
for T_probe in probe_horizons
haskey(results, T_probe) || continue
r = results[T_probe]
if r.status == "OK"
println(" T = $(T_probe) s: OK, $(r.n_sets) reach-sets, $(round(r.elapsed; digits=1))s wall")
else
println(" T = $(T_probe) s: FAILED")
end
end
# Save the longest-successful probe's envelope arrays for the app.
mat_out = joinpath(@__DIR__, "..", "..", "reachability", "reach_heatup_pj_result.mat")
saved = Dict{String, Any}()
saved["probe_horizons"] = collect(probe_horizons)
for T_probe in probe_horizons
haskey(results, T_probe) || continue
r = results[T_probe]
if r.status == "OK"
saved["T_$(Int(T_probe))_Tc_lo"] = r.Tc[1]
saved["T_$(Int(T_probe))_Tc_hi"] = r.Tc[2]
saved["T_$(Int(T_probe))_Tf_lo"] = r.Tf[1]
saved["T_$(Int(T_probe))_Tf_hi"] = r.Tf[2]
saved["T_$(Int(T_probe))_Tcold_lo"] = r.Tcold[1]
saved["T_$(Int(T_probe))_Tcold_hi"] = r.Tcold[2]
saved["T_$(Int(T_probe))_n_lo"] = r.n[1]
saved["T_$(Int(T_probe))_n_hi"] = r.n[2]
end
end
matwrite(mat_out, saved)
println("\nSaved envelope summaries to $mat_out")

113
code/scripts/validate_pj.jl Normal file
View File

@ -0,0 +1,113 @@
#!/usr/bin/env julia
#
# validate_pj.jl — quantify the prompt-jump approximation error.
#
# Two parallel sims of the same heatup scenario:
# (i) full 10-state PKE+T/H (pke_th_rhs!)
# (ii) reduced 9-state prompt-jump model (pke_th_rhs_pj!)
#
# After the prompt transient (~few hundred microseconds) the two
# trajectories should agree closely on T_c, T_f, T_cold, and on the
# reconstructed n vs the dynamic n. Quantify the error explicitly.
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using Printf
using LinearAlgebra
using OrdinaryDiffEq
using Plots
include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
include(joinpath(@__DIR__, "..", "src", "pke_th_rhs_pj.jl"))
include(joinpath(@__DIR__, "..", "controllers", "controllers.jl"))
plant = pke_params()
T_standby = plant.T_c0 - 33.333333
# Heatup scenario — same as sim_heatup.jl.
ref_heat = (; T_start=T_standby, T_target=plant.T_c0, ramp_rate=28/3600)
Q_sg = t -> 0.0
# Initial conditions — both at n=1e-3, same temperatures.
n0 = 1e-3
C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n0
x0_full = [n0; C0; T_standby; T_standby; T_standby]
x0_pj = [C0; T_standby; T_standby; T_standby]
tspan = (0.0, 3000.0)
# --- Full-state sim ---
function rhs_full!(dx, x, p, t)
u = ctrl_heatup(t, x, plant, ref_heat)
pke_th_rhs!(dx, x, t, plant, Q_sg, u)
end
prob_full = ODEProblem(rhs_full!, x0_full, tspan)
sol_full = solve(prob_full, Rodas5(); reltol=1e-8, abstol=1e-10)
# --- Prompt-jump sim ---
# ctrl_heatup expects 10-state x with T_f at index 8, T_c at 9.
# We pass an adapter that maps 9-state to the controller's expected layout.
function ctrl_heatup_pj(t, x_pj, plant_arg, ref_arg)
# Construct a fake 10-state with n placeholder; controller doesn't read n.
x10 = [0.0; x_pj[1:6]; x_pj[7]; x_pj[8]; x_pj[9]]
return ctrl_heatup(t, x10, plant_arg, ref_arg)
end
function rhs_pj!(dx, x, p, t)
u = ctrl_heatup_pj(t, x, plant, ref_heat)
pke_th_rhs_pj!(dx, x, t, plant, Q_sg, u)
end
prob_pj = ODEProblem(rhs_pj!, x0_pj, tspan)
sol_pj = solve(prob_pj, Rodas5(); reltol=1e-8, abstol=1e-10)
# --- Compare at sampled times ---
function get_n_full(sol, t) sol(t)[1] end
function get_T_c_full(sol, t) sol(t)[9] end
function get_T_f_full(sol, t) sol(t)[8] end
function get_T_cold_full(sol, t) sol(t)[10] end
function get_T_c_pj(sol, t) sol(t)[8] end
function get_T_f_pj(sol, t) sol(t)[7] end
function get_T_cold_pj(sol, t) sol(t)[9] end
function get_n_pj(sol, t)
x = sol(t)
u = ctrl_heatup_pj(t, x, plant, ref_heat)
pj_reconstruct_n(x, plant, u)
end
println("\n=== PJ vs full-state, heatup scenario ===")
println(" t [s] n_full n_pj |Δn|/n_full T_c err T_f err T_cold err")
for t_check in (1.0, 5.0, 10.0, 60.0, 300.0, 1200.0, 3000.0)
n_f = get_n_full(sol_full, t_check)
n_p = get_n_pj(sol_pj, t_check)
Tc_err = abs(get_T_c_full(sol_full, t_check) - get_T_c_pj(sol_pj, t_check))
Tf_err = abs(get_T_f_full(sol_full, t_check) - get_T_f_pj(sol_pj, t_check))
Tcold_err = abs(get_T_cold_full(sol_full, t_check) - get_T_cold_pj(sol_pj, t_check))
@printf " %6.1f %.3e %.3e %8.2e %.3e %.3e %.3e\n" t_check n_f n_p abs(n_f-n_p)/abs(n_f) Tc_err Tf_err Tcold_err
end
# --- Plot trajectory overlay ---
figdir = joinpath(@__DIR__, "..", "..", "docs", "figures")
isdir(figdir) || mkpath(figdir)
CtoF(T) = T*9/5 + 32
t_grid = range(0, 3000, length=600)
n_full_arr = [get_n_full(sol_full, t) for t in t_grid]
n_pj_arr = [get_n_pj(sol_pj, t) for t in t_grid]
Tc_full_arr = [get_T_c_full(sol_full, t) for t in t_grid]
Tc_pj_arr = [get_T_c_pj(sol_pj, t) for t in t_grid]
p_n = plot(t_grid ./ 60, n_full_arr, lw=2, color=:blue, label="full-state",
xlabel="t [min]", ylabel="n", title="Power: full vs PJ")
plot!(p_n, t_grid ./ 60, n_pj_arr, lw=1.5, ls=:dash, color=:red, label="prompt-jump")
p_Tc = plot(t_grid ./ 60, CtoF.(Tc_full_arr), lw=2, color=:blue, label="full-state",
xlabel="t [min]", ylabel="T_avg [F]", title="T_avg: full vs PJ")
plot!(p_Tc, t_grid ./ 60, CtoF.(Tc_pj_arr), lw=1.5, ls=:dash, color=:red, label="prompt-jump")
fig = plot(p_n, p_Tc, layout=(1,2), size=(1200, 450),
plot_title="PJ vs full-state, heatup scenario (3000 s)")
savefig(fig, joinpath(figdir, "validate_pj_heatup.png"))
println("\nFigure: $figdir/validate_pj_heatup.png")

104
code/src/pke_th_rhs_pj.jl Normal file
View File

@ -0,0 +1,104 @@
"""
pke_th_rhs_pj!(dx, x, t, plant, Q_sg, u)
Prompt-jump (singular-perturbation) reduced-order PKE + thermal-hydraulics.
The full PKE has prompt-neutron timescale Λ ≈ 10⁻⁴ s, while thermal
dynamics evolve on 1–100 s. For reach analysis on horizons longer
than ~10 s, the prompt transient is irrelevant to safety but its
~50 µs timescale forces TMJets to take ~ms steps, exhausting the
step budget.
Singular-perturbation reduction: set dn/dt ≈ 0 in the neutron-balance
equation:
0 = (ρ - β)/Λ · n + Σ λᵢ Cᵢ
→ n ≈ Λ Σ λᵢ Cᵢ / (β - ρ)
Valid when β - ρ > 0 (sub-prompt-critical), which holds for normal
operation where ρ varies on the scale of pcm vs β = 650 pcm.
State drops from 10 to 9:
x = [C₁..C₆, T_f, T_c, T_cold] (no n; n is computed algebraically)
The C dynamics still depend on n; we substitute the algebraic form,
giving:
dCᵢ/dt = (βᵢ/Λ) · n - λᵢ · Cᵢ
= βᵢ · Σⱼ λⱼ Cⱼ / (β - ρ) - λᵢ · Cᵢ
This introduces a rational-function nonlinearity (1 / (β - ρ)) into
both the precursor balance and the fuel-temperature equation. Taylor
models handle it as long as the denominator stays bounded away from
zero — a non-issue under normal operation.
Reach-analysis benefit: removes the Λ⁻¹ stiffness, so step size is
limited only by the slowest-precursor + thermal timescales.
Soundness cost: the prompt transient (∼50 µs after a reactivity
perturbation) is no longer captured. For safety claims on hours-long
horizons this is acceptable; for prompt-supercritical excursions it is
fundamentally wrong (the algebraic form blows up as ρ → β).
"""
function pke_th_rhs_pj!(dx, x, t, plant, Q_sg, u)
# Unpack 9-state x = [C1..C6, T_f, T_c, T_cold]
C1, C2, C3, C4, C5, C6 = x[1], x[2], x[3], x[4], x[5], x[6]
T_f = x[7]
T_c = x[8]
T_cold = x[9]
T_hot = 2 * T_c - T_cold
# Total reactivity (controller u plus T-feedback).
rho = u + plant.alpha_f * (T_f - plant.T_f0) +
plant.alpha_c * (T_c - plant.T_c0)
# Prompt-jump algebraic n.
sum_lam_C = plant.lambda_i[1]*C1 + plant.lambda_i[2]*C2 +
plant.lambda_i[3]*C3 + plant.lambda_i[4]*C4 +
plant.lambda_i[5]*C5 + plant.lambda_i[6]*C6
denom = plant.beta - rho
n = plant.Lambda * sum_lam_C / denom
inv_factor = sum_lam_C / denom # used by precursor balance
# Precursor balance under PJ.
dx[1] = plant.beta_i[1] * inv_factor - plant.lambda_i[1] * C1
dx[2] = plant.beta_i[2] * inv_factor - plant.lambda_i[2] * C2
dx[3] = plant.beta_i[3] * inv_factor - plant.lambda_i[3] * C3
dx[4] = plant.beta_i[4] * inv_factor - plant.lambda_i[4] * C4
dx[5] = plant.beta_i[5] * inv_factor - plant.lambda_i[5] * C5
dx[6] = plant.beta_i[6] * inv_factor - plant.lambda_i[6] * C6
# Thermal-hydraulics (n appears algebraically in the fuel eq.).
dx[7] = (plant.P0 * n - plant.hA * (T_f - T_c)) / (plant.M_f * plant.c_f)
dx[8] = (plant.hA * (T_f - T_c) - 2 * plant.W * plant.c_c * (T_c - T_cold)) /
(plant.M_c * plant.c_c)
dx[9] = (plant.W * plant.c_c * (T_hot - T_cold) - Q_sg(t)) /
(plant.M_sg * plant.c_c)
return nothing
end
"""
pke_initial_conditions_pj(plant; n_target=1.0)
Reduced-state initial condition for the prompt-jump model. Returns
the 9-element vector [C1..C6, T_f, T_c, T_cold] consistent with
n = `n_target` at full-power equilibrium.
"""
function pke_initial_conditions_pj(plant; n_target=1.0)
C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n_target
return [C0; plant.T_f0; plant.T_c0; plant.T_cold0]
end
"""
pj_reconstruct_n(x, plant, u)
Algebraic n from the PJ reduced state. For diagnostics / plotting.
"""
function pj_reconstruct_n(x, plant, u)
T_f = x[7]
T_c = x[8]
rho = u + plant.alpha_f * (T_f - plant.T_f0) +
plant.alpha_c * (T_c - plant.T_c0)
sum_lam_C = sum(plant.lambda_i .* x[1:6])
return plant.Lambda * sum_lam_C / (plant.beta - rho)
end

BIN
docs/figures/validate_pj_heatup.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 59 KiB

218
journal/entries/2026-04-20-overnight-prompt-jump.tex Normal file
View File

@ -0,0 +1,218 @@
% ---------------------------------------------------------------------------
% 2026-04-20 — overnight session: prompt-jump reduction + nonlinear reach
% A-style invention-log entry, written in real-time during the session.
% ---------------------------------------------------------------------------
\session{2026-04-20 (overnight)}{open-ended autonomous session}{Implement
the singular-perturbation (prompt-jump) reduction of the PKE+T/H model.
Validate it against the full 10-state. Re-run TMJets nonlinear reach
on heatup and find the new horizon wall. Extend the Pluto app to read
reach results live. Document everything for review in the morning.}
\section{2026-04-20 (overnight) — Prompt-jump nonlinear reach}
\label{sec:20260420-overnight}
\subsection*{Origin}
The 2026-04-20 evening session ended with TMJets working on the full
10-state heatup at $T = 10$~\unit{\second} but exhausting its 50{,}000-step
budget by $T = 60$~\unit{\second}. Diagnosis: prompt-neutron timescale
$\Lambda = 10^{-4}$~\unit{\second} forces $\sim$1~\unit{\milli\second}
adaptive steps to bound Taylor remainder. Over hours, infeasible.
The known remedy: \emph{singular-perturbation reduction} — set
$\dot n = 0$ and solve algebraically for $n$, removing the prompt
timescale from the dynamic state. Standard reactor-kinetics move,
documented in textbooks (Hetrick \emph{Dynamics of Nuclear Reactors},
ch.\ 4; Ott \& Neuhold). Auto mode active; Dane's instruction at
session start: ``take a big fat overnight rip as far as you can on the
prompt jump assumption and doing the reachability and app buildout.
Document things in the journal and we'll review in the morning.''
\subsection*{Part 1: The prompt-jump derivation}
\begin{derivation}
Starting from the 10-state PKE+T/H system, focus on the neutron-balance
equation:
$$\dot n = \frac{\rho - \beta}{\Lambda} n + \sum_{i=1}^{6} \lambda_i C_i.$$
The prompt-neutron generation time $\Lambda \sim 10^{-4}$~\unit{\second}
makes the first term \emph{very fast} relative to the precursor and
thermal dynamics (precursor timescales 0.3 to 80~\unit{\second}; thermal
$\sim$10--100~\unit{\second}). A standard regular-perturbation argument
(Hetrick, ch.~4) shows that on timescales $\gg \Lambda$, the prompt
term equilibrates rapidly and we can set
$$\dot n \approx 0 \quad \Longrightarrow \quad
\frac{\rho - \beta}{\Lambda} n + \sum_i \lambda_i C_i = 0.$$
Solving for $n$:
$$\boxed{\;n_{\mathrm{PJ}}(C, \rho) = \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho}\;}$$
valid when $\beta - \rho > 0$, i.e.\ sub-prompt-critical. For our
heatup controller $\rho = K_p \cdot e$ with $K_p e \ll \beta$, so the
denominator is well bounded away from zero.
Substituting back into the precursor and fuel equations:
\begin{align*}
\dot C_i &= \frac{\beta_i}{\Lambda} n_{\mathrm{PJ}} - \lambda_i C_i
= \frac{\beta_i \sum_j \lambda_j C_j}{\beta - \rho} - \lambda_i C_i \\
\dot T_f &= \frac{P_0 \, n_{\mathrm{PJ}} - hA(T_f - T_c)}{M_f c_f}
= \frac{P_0 \Lambda \sum_j \lambda_j C_j / (\beta - \rho) - hA(T_f - T_c)}{M_f c_f}.
\end{align*}
The state vector drops from 10 to 9: $x = [C_1, \ldots, C_6, T_f, T_c, T_{\mathrm{cold}}]^\top$.
The dynamics gain a rational nonlinearity ($1/(\beta - \rho)$). The
fastest dynamic timescale becomes $1/\lambda_6 = 0.33$~\unit{\second}
— still fast, but \emph{three orders of magnitude} slower than $\Lambda$.
\textbf{Soundness cost:} the prompt transient (the $\sim$50~\unit{\micro\second}
adjustment of $n$ after a step in $\rho$) is no longer captured. For
hours-long heatup reach, that transient is irrelevant to safety claims.
For prompt-supercritical regimes ($\rho \to \beta$) the algebraic
formula diverges and the reduction is invalid — but those regimes are
themselves accident-class, outside the scope of normal-operation reach.
\end{derivation}
\subsection*{Part 2: Implementation}
Two new files in \texttt{code/}:
\begin{itemize}
\item \texttt{src/pke\_th\_rhs\_pj.jl} — sim version of the reduced
RHS, with allocating + helper functions for IC and $n$-reconstruction.
\item \texttt{scripts/validate\_pj.jl} — side-by-side sim of full
vs.\ reduced PKE on the heatup scenario.
\end{itemize}
The reduced RHS is structurally identical to the full one with two
differences: (a) no $\dot n$ equation; (b) $n$ inside the precursor and
fuel-temperature equations is replaced by $n_{\mathrm{PJ}}(C, \rho)$,
introducing the rational denominator.
\subsection*{Part 3: Validation against full-state}
\texttt{validate\_pj.jl} runs both models from the same heatup IC
($n_0 = 10^{-3}$, $T = T_{\mathrm{standby}}$ everywhere) for 50 minutes
and tabulates pointwise error.
\begin{lstlisting}[style=terminal]
=== PJ vs full-state, heatup scenario ===
t [s] n_full n_pj |Δn|/n_full T_c err T_f err T_cold err
1.0 1.000e-03 1.000e-03 8.32e-07 4.839e-09 1.718e-08 6.642e-10
5.0 1.000e-03 1.000e-03 3.08e-06 3.970e-08 9.392e-08 1.921e-08
10.0 1.001e-03 1.001e-03 5.59e-06 1.295e-07 2.320e-07 7.945e-08
60.0 1.017e-03 1.018e-03 3.70e-05 3.826e-06 4.534e-06 3.446e-06
300.0 1.310e-03 1.311e-03 3.77e-04 1.867e-04 1.960e-04 1.816e-04
1200.0 3.414e-03 3.410e-03 1.02e-03 2.177e-03 2.111e-03 2.213e-03
3000.0 3.248e-03 3.250e-03 5.03e-04 7.166e-03 7.197e-03 7.149e-03
\end{lstlisting}
\textbf{Maximum relative error on $n$ over 3000~\unit{\second}: 0.10\%}
(at $t = 1200$~\unit{\second}). Maximum temperature error: 7~\unit{\milli\kelvin}.
The PJ approximation is excellent — the absolute errors are far below
any physical safety margin.
The PJ trajectory is essentially indistinguishable from full-state on
the heatup timescale (\cref{fig:validate-pj}).
\begin{figure}[h]
\centering
\includegraphics[width=0.95\linewidth]{validate_pj_heatup.png}
\caption{Full-state (blue) vs.\ prompt-jump (red dashed) sims of the
same heatup scenario. Power $n$ (left) and $T_{\mathrm{avg}}$
(right) overlay almost perfectly across 50~\unit{\minute}. The
difference is invisible at this scale — peak relative error on $n$
is 0.1\%. This is the empirical evidence that the singular-perturbation
reduction is sound for this class of slow heatup transients.}
\label{fig:validate-pj}
\end{figure}
\subsection*{Part 4: Nonlinear reach with the PJ model}
\apass{Results are populating as TMJets runs in the background. Final
horizon and step counts will be filled in here once the longest probe
returns. Initial expectation: 60~\unit{\second} should pass trivially;
1800~\unit{\second} is the real test; 5400 / 18000 the dream.}
The PJ reach script is \texttt{code/scripts/reach\_heatup\_pj.jl}.
Same Taylor-model machinery (\texttt{TMJets}, \texttt{@taylorize},
augmented time state) as the failed full-state version, but the RHS
operates on the 9-state PJ system (10D with augmented time) and
includes the rational $1/(\beta - \rho)$ in two places. Probe
horizons: 60, 300, 1800, 5400~\unit{\second}.
\begin{decision}
TMJets settings: \texttt{orderT=4}, \texttt{orderQ=2}, \texttt{abstol=1e-9},
\texttt{maxsteps=100\,000}. \texttt{abstol} is one order looser than
the full-state attempt — the PJ RHS has a rational nonlinearity that
narrows the Taylor remainder convergence radius slightly, and we don't
need 1e-10 precision for envelope tracking on a tube that's already
several Kelvin wide.
\end{decision}
\apass{Results section will be populated below as probes complete.
Currently TMJets is precompiling (~5--15 minutes for the new
\texttt{@taylorize}'d RHS).}
\subsection*{Part 5: App buildout}
While the reach is running, extended the Pluto predicate explorer
with three new sections:
\begin{itemize}
\item \textbf{Live reach-result ingestion} (§9b): reads
\texttt{reachability/reach\_operation\_result.mat} (saved by
\texttt{reach\_operation.jl}) and renders per-halfspace margins
live, replacing the hand-maintained traceability table.
\item \textbf{2D projection chooser} (§9c): pick any two state
coordinates from $\{n, T_f, T_c, T_{\mathrm{cold}}\}$ and see
the operating polytope with the reach-tube envelope as a red
rectangle overlay.
\item \textbf{PJ heatup reach overlay} (§9d): if \texttt{reach\_heatup\_pj\_result.mat}
exists, display the envelope summary.
\end{itemize}
Added \texttt{MAT.jl} to the app's \texttt{Project.toml}. Read-only
v1 still — sliders preview UX without writing back.
\subsection*{Soundness ledger update}
\begin{decision}
The PJ reduction shifts the soundness story:
\textbf{Before:} linear reach was a sound over-approximation of the
linearized closed-loop, but the linearization was an unbounded
approximation of the nonlinear plant. Net: \emph{approximate, not
sound} for the plant.
\textbf{After:} TMJets nonlinear reach with PJ is a sound
over-approximation of the \emph{prompt-jump-reduced} nonlinear plant.
The PJ reduction itself introduces a controlled approximation
(0.1\% error on $n$, mK on $T$, validated empirically over 50
minutes). Net: \emph{$\epsilon_{\mathrm{PJ}}$-approximate but otherwise
sound}, where $\epsilon_{\mathrm{PJ}}$ is bounded.
This is qualitatively better. The remaining gap (PJ approximation
error) can be characterized by the validation experiment, which we have.
The next step toward full soundness would be a Tikhonov-style
singular-perturbation theorem application giving a closed-form
$\mathcal{O}(\Lambda)$ error bound, but the empirical bound is
defensible for the prelim demo.
\end{decision}
\subsection*{Open at close}
\apass{This entry is being written in parallel with the running
reach. Final results to be filled in below as TMJets returns. If
TMJets completes the 5-hour horizon, the heatup reach-avoid obligation
is discharged (modulo PJ + saturation caveats). If it stops earlier,
identify the new wall and propose the next reduction.}
\begin{itemize}
\item Polytopic / SOS barriers — still the only path to a tight
analytic certificate; quadratic Lyapunov is structurally
defeated regardless of model order.
\item Saturation as explicit hybrid sub-mode — still pending,
independent of PJ.
\item Parametric $\alpha$ uncertainty — still pending.
\item Tikhonov / regular-perturbation $\mathcal{O}(\Lambda)$ error
bound on PJ.
\item Per-mode reach for shutdown and scram (now feasible with PJ).
\end{itemize}

2
journal/journal.tex
View File

@ -72,5 +72,7 @@ Each limitation ties to a plan or an open question.
\input{entries/2026-04-20-predicates-boundaries-julia-nonlinear.tex}
\newpage
\input{entries/2026-04-20-evening-mega-session.tex}
\newpage
\input{entries/2026-04-20-overnight-prompt-jump.tex}
\end{document}