julia: nonlinear heatup reach — 10s horizon works, longer fails on stiffness

First working nonlinear reach artifact for the PWR model. TMJets
Taylor-model scheme on the full 10-state closed-loop (unsaturated
ctrl_heatup, ramp reference via augmented time state x[11]).

Status:
  T=10s    : SUCCESS. 10583 reach-sets. T_c envelope [274.45, 295] C,
             n envelope [-5.2e-4, 5.01e-3]. Over-approximation visible
             (n can't be negative physically) but tube is sound and
             bounded.
  T=60s+   : FAILED. Exhausts 50k step budget then hits NaN in
             precursor-decay term.

Root cause: prompt-neutron stiffness. Lambda=1e-4s forces TMJets'
adaptive stepper to ~1ms steps to resolve fast dynamics. 10583 steps
for 10s of sim time means we get ~10s/50000 = 2s horizon max before
step budget exhausts — inadequate for heatup's 5-hour obligation.

Remedy (next session): singular-perturbation reduction of the
neutronics. Treat n as quasi-static algebraic function of (T, C, rho)
rather than a dynamic state. Replaces stiff dn/dt with algebraic
relation, removes fast timescale from reach problem. Standard in
reactor-kinetics reach literature.

What this does prove:
  - Julia/TMJets framework works for this system (previous
    scaling-issue failure is gone with @taylorize'd RHS).
  - Bilinear n*rho term handled correctly by Taylor models.
  - Ramp reference via augmented time state x[11] is a workable
    pattern for time-varying controller references in reach.

What this does NOT prove:
  - Anything about heatup safety — 10s horizon is nowhere near the
    mode's 5-hour obligation.

Includes sim_heatup.jl, a Rodas5 baseline using the same @taylorize-
able RHS form, for cross-validation of the reach tube against a
nominal trajectory once longer horizons are reachable.

WALKTHROUGH.md updated with the finding.

Hacker-Split: got partway up the Julia reach ladder, identified the
physical bottleneck (stiffness), named the fix (reduced-order PKE).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

julia-port/Project.toml

@@ -1,6 +1,7 @@
 authors = ["Dane Sabo <yourstruly@danesabo.com>"]
 
 [deps]
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"

julia-port/controllers/ctrl_heatup_unsat.jl

@@ -0,0 +1,20 @@
+"""
+    ctrl_heatup_unsat(t, x, plant, ref)
+
+Heatup controller WITHOUT saturation.  Identical to `ctrl_heatup` but the
+`clamp()` is removed — u can be anything the feedback-linearization +
+ramped-reference P wants.
+
+Used as the starting point for nonlinear reach analysis.  The full
+`ctrl_heatup` adds saturation on top of this, which turns the controller
+into a 3-mode piecewise-affine system and is hybrid-reach territory.
+"""
+function ctrl_heatup_unsat(t, x, plant, ref)
+    Kp = 1e-4
+    T_f   = x[8]
+    T_avg = x[9]
+    u_ff = -plant.alpha_f * (T_f   - plant.T_f0) -
+            plant.alpha_c * (T_avg - plant.T_c0)
+    T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target)
+    return u_ff + Kp * (T_ref - T_avg)
+end

julia-port/scripts/reach_heatup_nonlinear.jl

@@ -0,0 +1,176 @@
+#!/usr/bin/env julia
+#
+# reach_heatup_nonlinear.jl — nonlinear reach on heatup mode via TMJets.
+#
+# STATUS: partial success. 10 s horizon works; longer horizons fail.
+#
+# Full 10-state nonlinear reach of the PKE + thermal-hydraulics
+# closed-loop is limited to ~10 s horizons by prompt-neutron stiffness
+# (Lambda = 1e-4 s).  TMJets' adaptive step size shrinks to ~1 ms to
+# resolve prompt dynamics, then hits numerical precision floor and
+# produces NaN around t ~ 10 s.  10,583 steps for 10 s of sim time.
+#
+# For heatup's 5-hour horizon (18 million prompt timescales), we need
+# singular-perturbation reduction of the neutronics: treat n as a
+# quasi-static algebraic function of (T, C, rho) rather than a dynamic
+# state.  The reduced-order PKE replaces the stiff dn/dt equation with
+# an algebraic relation, removing the fast timescale from the reach
+# problem.
+#
+# The 10 s result is a validation of the Taylor-model approach, not
+# a usable safety proof.  The machinery works; the model needs
+# reduction before the horizon is meaningful.
+#
+# Saturation disabled (ctrl_heatup_unsat structure), all constants
+# inlined so @taylorize can generate a specialized Taylor-model RHS.
+# Time carried as augmented state x[11].
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+
+using LinearAlgebra
+using ReachabilityAnalysis, LazySets
+using JSON
+
+# --- Plant constants (inlined; must match pke_params.m) ---
+const LAMBDA  = 1e-4
+const BETA_1  = 0.000215
+const BETA_2  = 0.001424
+const BETA_3  = 0.001274
+const BETA_4  = 0.002568
+const BETA_5  = 0.000748
+const BETA_6  = 0.000273
+const BETA    = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6
+const LAM_1   = 0.0124
+const LAM_2   = 0.0305
+const LAM_3   = 0.111
+const LAM_4   = 0.301
+const LAM_5   = 1.14
+const LAM_6   = 3.01
+
+const P0   = 1e9
+const M_F  = 50000.0
+const C_F  = 300.0
+const M_C  = 20000.0
+const C_C  = 5450.0
+const HA   = 5e7
+const W_M  = 5000.0
+const M_SG = 30000.0
+
+const ALPHA_F = -2.5e-5
+const ALPHA_C = -1e-4
+
+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
+
+# --- Controller constants ---
+const T_STANDBY    = T_C0 - 33.333333
+const RAMP_RATE_CS = 28.0 / 3600
+const KP_HEATUP    = 1e-4
+
+# --- Taylorized RHS with augmented time (x[11] = t, dx[11] = 1) ---
+# TMJets handles augmented time better than accessing the solver's t arg
+# directly inside the Taylor-model rewriting.
+@taylorize function rhs_heatup_taylor!(dx, x, p, t)
+    # Inlined cancellation: rho_total = Kp * (T_ref - T_c) after cancellation.
+    # T_ref uses x[11] as the time-state (no min; valid while T_ref < T_c0).
+    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[11] - x[9])
+
+    # Neutronics (explicit — no loops, no function calls)
+    dx[1] = (rho - BETA) / LAMBDA * x[1] +
+            LAM_1*x[2] + LAM_2*x[3] + LAM_3*x[4] +
+            LAM_4*x[5] + LAM_5*x[6] + LAM_6*x[7]
+    dx[2] = (BETA_1 / LAMBDA) * x[1] - LAM_1 * x[2]
+    dx[3] = (BETA_2 / LAMBDA) * x[1] - LAM_2 * x[3]
+    dx[4] = (BETA_3 / LAMBDA) * x[1] - LAM_3 * x[4]
+    dx[5] = (BETA_4 / LAMBDA) * x[1] - LAM_4 * x[5]
+    dx[6] = (BETA_5 / LAMBDA) * x[1] - LAM_5 * x[6]
+    dx[7] = (BETA_6 / LAMBDA) * x[1] - LAM_6 * x[7]
+
+    # Thermal-hydraulics
+    dx[8]  = (P0 * x[1] - HA * (x[8] - x[9])) / (M_F * C_F)
+    dx[9]  = (HA * (x[8] - x[9]) - 2 * W_M * C_C * (x[9] - x[10])) / (M_C * C_C)
+    dx[10] = (2 * W_M * C_C * (x[9] - x[10])) / (M_SG * C_C)
+
+    # Time (augmented state for reach; dt/dt = 1)
+    dx[11] = one(x[1])
+    return nothing
+end
+
+# --- Build X_entry 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
+
+x_lo = [n_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]   # augmented time
+
+x_hi = [n_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 (saturation disabled, @taylorize RHS) ===")
+println("  X_entry: n ∈ [$(n_lo), $(n_hi)], T_c ∈ [$(T_c_lo), $(T_c_hi)] C")
+
+for T_probe in (10.0, 60.0, 300.0)
+    println("\n--- Probe T = $T_probe s ---")
+    sys  = BlackBoxContinuousSystem(rhs_heatup_taylor!, 11)
+    prob = InitialValueProblem(sys, X0)
+    try
+        alg = TMJets(orderT=4, orderQ=2, abstol=1e-10, maxsteps=50000)
+        sol = solve(prob; T=T_probe, alg=alg)
+        flow = flowpipe(sol)
+        n_sets = length(flow)
+        println("  TMJets: $n_sets reach-sets")
+
+        # Overapproximate each Taylor-model reach-set as a hyperrectangle
+        # so we can query envelopes.  This is standard LazySets usage.
+        flow_hr = overapproximate(flow, Hyperrectangle)
+        n_sets_hr = length(flow_hr)
+
+        # Envelope per state over the whole flowpipe.
+        n_lo_env  = minimum(low(set(R), 1)  for R in flow_hr)
+        n_hi_env  = maximum(high(set(R), 1) for R in flow_hr)
+        Tc_lo_env = minimum(low(set(R), 9)  for R in flow_hr)
+        Tc_hi_env = maximum(high(set(R), 9) for R in flow_hr)
+        Tf_lo_env = minimum(low(set(R), 8)  for R in flow_hr)
+        Tf_hi_env = maximum(high(set(R), 8) for R in flow_hr)
+        Tcold_lo  = minimum(low(set(R), 10) for R in flow_hr)
+        Tcold_hi  = maximum(high(set(R), 10) for R in flow_hr)
+        t_final   = maximum(high(set(R), 11) for R in flow_hr)
+
+        println("  t reached:    $(round(t_final; digits=1)) s of $T_probe")
+        println("  n envelope:   [$(round(n_lo_env; sigdigits=4)), $(round(n_hi_env; 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 env:   [$(round(Tcold_lo; digits=2)), $(round(Tcold_hi; digits=2))] C")
+    catch err
+        msg = sprint(showerror, err)
+        println("  FAILED: ", first(msg, 300))
+    end
+end

julia-port/scripts/sim_heatup.jl

@@ -0,0 +1,73 @@
+#!/usr/bin/env julia
+#
+# sim_heatup.jl — nominal heatup trajectory sim in Julia.
+#
+# Closed-loop nonlinear simulation using the SAME RHS form as
+# reach_heatup_nonlinear.jl (augmented time x[11], inlined constants)
+# but via OrdinaryDiffEq, not TMJets.  Purpose: verify the @taylorize-able
+# RHS matches the MATLAB heatup baseline.  Also generates a nominal
+# trajectory that reach tubes can be checked against.
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+
+using OrdinaryDiffEq
+
+# Constants (match reach_heatup_nonlinear.jl).
+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
+const ALPHA_F, ALPHA_C = -2.5e-5, -1e-4
+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
+
+function rhs_heatup_sim!(dx, x, p, t)
+    rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[11] - x[9])
+    dx[1] = (rho - BETA) / LAMBDA * x[1] +
+            LAM_1*x[2] + LAM_2*x[3] + LAM_3*x[4] +
+            LAM_4*x[5] + LAM_5*x[6] + LAM_6*x[7]
+    dx[2] = (BETA_1 / LAMBDA) * x[1] - LAM_1 * x[2]
+    dx[3] = (BETA_2 / LAMBDA) * x[1] - LAM_2 * x[3]
+    dx[4] = (BETA_3 / LAMBDA) * x[1] - LAM_3 * x[4]
+    dx[5] = (BETA_4 / LAMBDA) * x[1] - LAM_4 * x[5]
+    dx[6] = (BETA_5 / LAMBDA) * x[1] - LAM_5 * x[6]
+    dx[7] = (BETA_6 / LAMBDA) * x[1] - LAM_6 * x[7]
+    dx[8]  = (P0 * x[1] - HA * (x[8] - x[9])) / (M_F * C_F)
+    dx[9]  = (HA * (x[8] - x[9]) - 2 * W_M * C_C * (x[9] - x[10])) / (M_C * C_C)
+    dx[10] = (2 * W_M * C_C * (x[9] - x[10])) / (M_SG * C_C)
+    dx[11] = 1.0
+    return nothing
+end
+
+# IC = midpoint of X_entry polytope.
+n0 = 0.001
+C0 = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA), BETA_3/(LAM_3*LAMBDA),
+      BETA_4/(LAM_4*LAMBDA), BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n0
+x0 = [n0; C0; T_STANDBY + 5; T_STANDBY + 6; T_STANDBY + 4; 0.0]
+
+tspan = (0.0, 300.0)
+
+# Rodas5 is stiff-aware.
+prob = ODEProblem(rhs_heatup_sim!, x0, tspan)
+sol  = solve(prob, Rodas5(); reltol=1e-8, abstol=1e-10)
+
+println("\n=== Heatup nominal trajectory (saturation disabled, Julia) ===")
+println("  t=0       : n=$(round(sol.u[1][1]; sigdigits=3))  T_c=$(round(sol.u[1][9]; digits=2)) C")
+for t_check in (10.0, 60.0, 120.0, 300.0)
+    u = sol(t_check)
+    println("  t=$(lpad(Int(t_check), 4))s    : n=$(round(u[1]; sigdigits=3))  T_c=$(round(u[9]; digits=2)) C  T_f=$(round(u[8]; digits=2)) C")
+end

reachability/WALKTHROUGH.md

@@ -448,11 +448,22 @@ vibes, not a proof.
 
 ### What's next
 
-1. **Nonlinear reach on heatup via JuliaReach.** First pass: saturation
+1. **Nonlinear reach on heatup via JuliaReach — partial result in.**
-   disabled, controller free to apply any u. See how Taylor models
+   `../julia-port/scripts/reach_heatup_nonlinear.jl` ran successfully
-   handle the bilinear `n·e` term. If TMJets produces a tight tube,
+   to T=10 s on the full 10-state nonlinear closed-loop (saturation
-   we have a soundness story for transition modes. (The Julia port
+   disabled). Got 10,583 reach-sets with envelope T_c ∈ [274.45, 295]
-   scaffolding is in `../julia-port/`.)
+   and n ∈ [-5e-4, 5e-3]. Nonlinear reach is *functional* in Julia —
+   the framework, `@taylorize` RHS, and bilinearity handling all work.
+   **But the horizon is limited to ~10 s by prompt-neutron stiffness:**
+   TMJets' adaptive stepper shrinks to ~1 ms per step to resolve
+   Λ = 10⁻⁴ s prompt-neutron dynamics, then hits numerical precision
+   floor and produces NaN around t = 10-20 s. For heatup's 5-hour
+   horizon (18 million prompt timescales) we need a **reduced-order
+   model**: singular-perturbation elimination of the fast neutronics
+   (treat n quasi-statically as an algebraic function of thermal
+   states + precursors). Standard in reactor-kinetics reach literature.
+   This is the actual next step before nonlinear reach can produce a
+   usable safety proof for heatup.
 2. **Once nonlinear reach works: add saturation as hybrid sub-mode.**
 3. **Parametric α reach** once the framework can handle it.
 4. **Shutdown and scram reach** (trivial cases, should be quick once