diff --git a/app/Project.toml b/app/Project.toml index 3b26fbb..6d14573 100644 --- a/app/Project.toml +++ b/app/Project.toml @@ -2,6 +2,7 @@ authors = ["Dane Sabo "] [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" diff --git a/app/predicate_explorer.jl b/app/predicate_explorer.jl index 5ed5371..d1dacc8 100644 --- a/app/predicate_explorer.jl +++ b/app/predicate_explorer.jl @@ -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 diff --git a/code/scripts/reach_heatup_pj.jl b/code/scripts/reach_heatup_pj.jl new file mode 100644 index 0000000..bf49613 --- /dev/null +++ b/code/scripts/reach_heatup_pj.jl @@ -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") diff --git a/code/scripts/validate_pj.jl b/code/scripts/validate_pj.jl new file mode 100644 index 0000000..894af04 --- /dev/null +++ b/code/scripts/validate_pj.jl @@ -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") diff --git a/code/src/pke_th_rhs_pj.jl b/code/src/pke_th_rhs_pj.jl new file mode 100644 index 0000000..1df66e8 --- /dev/null +++ b/code/src/pke_th_rhs_pj.jl @@ -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 diff --git a/docs/figures/validate_pj_heatup.png b/docs/figures/validate_pj_heatup.png new file mode 100644 index 0000000..d5c6459 Binary files /dev/null and b/docs/figures/validate_pj_heatup.png differ diff --git a/journal/entries/2026-04-20-overnight-prompt-jump.tex b/journal/entries/2026-04-20-overnight-prompt-jump.tex new file mode 100644 index 0000000..05a4ac2 --- /dev/null +++ b/journal/entries/2026-04-20-overnight-prompt-jump.tex @@ -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} diff --git a/journal/journal.tex b/journal/journal.tex index 6a19ba0..bdb0c36 100644 --- a/journal/journal.tex +++ b/journal/journal.tex @@ -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}