PWR-HYBRID-3/journal/entries/2026-04-20-predicates-boundaries-julia-nonlinear.tex

% ---------------------------------------------------------------------------
% 2026-04-20 — Predicates restructure, mode taxonomy, Julia nonlinear reach
% Deep / A-style invention-log entry.
% ---------------------------------------------------------------------------

\session{2026-04-20}{afternoon session, $\sim$3 hr}{Walk through the
2026-04-17 work with a critical eye, surface three structural errors
hiding in it, restructure the predicates JSON to fix them, formalize
the per-mode reach-obligation taxonomy (equilibrium vs.\ transition),
and take the first swing at Julia nonlinear reach.  The Julia reach
framework works but is limited to $\sim$10-second horizons by
prompt-neutron stiffness.  The remedy (singular-perturbation reduction)
is identified and deferred.}

\section{2026-04-20 — Predicates restructure, mode taxonomy, nonlinear reach}
\label{sec:20260420-afternoon}

\subsection*{How this session started}

Walking through the 2026-04-17 work slide-by-slide with a critical
reviewer looking for inconsistencies.  Each section of the linear
reach story had a latent bug the first pass missed; by the end of the
walkthrough three structural errors were exposed and fixed.

\subsection*{Error 1: the barrier was checking the wrong invariant}

The 2026-04-17 barrier code compared the Lyapunov-ellipsoid reach
against a \emph{symmetric} slab $|T_c - T_{c0}| \leq 2.78$~\unit{\celsius}
— the operational deadband \texttt{t\_avg\_in\_range}.  But that
predicate is used by the DRC for \emph{mode transitions}: crossing it
triggers heatup$\to$operation, not damage.  The barrier should be
checking the \emph{safety limits} — one-sided halfspaces reflecting
actual trip setpoints and physical damage thresholds.  These are not
symmetric:
\begin{itemize}
  \item $T_f \leq 1200$~\unit{\celsius} (fuel centerline design limit).
  \item $T_c \leq 320$~\unit{\celsius} (high-$T_{\mathrm{avg}}$ trip).
  \item $T_c \geq 280$~\unit{\celsius} (low-$T_{\mathrm{avg}}$ trip).
  \item $n \leq 1.15$ (high-flux trip).
  \item $n \geq 0.15$ (operation range only).
  \item $T_{\mathrm{cold}} \leq T_{\mathrm{cold0}} + 15$ (subcooling).
\end{itemize}

\begin{decision}
Restructure \texttt{reachability/predicates.json} into three categories:
\begin{description}
  \item[\texttt{operational\_deadbands}] soft bands used by DRC for
        mode switching.  \texttt{t\_avg\_above\_min},
        \texttt{t\_avg\_in\_range}, \texttt{p\_above\_crit}.
  \item[\texttt{safety\_limits}] hard one-sided halfspaces.  Six rows
        as listed above, plus new heatup-rate ones (see Error 3).
  \item[\texttt{mode\_invariants}] \texttt{inv1\_holds},
        \texttt{inv2\_holds}: declared as \emph{conjunctions} of
        named safety limits, so changing a limit propagates.
\end{description}
Safety proofs target \texttt{mode\_invariants}; DRC transitions
target \texttt{operational\_deadbands}.  The distinction is
load-bearing.
\end{decision}

\subsection*{Error 2: the "2364$\times$ bound" should have included LQR, and didn't tell us}

After the restructure, re-running \texttt{barrier\_lyapunov.m} against
\texttt{inv2\_holds} produced:

\begin{lstlisting}[style=terminal]
=== Lyapunov barrier vs inv2_holds halfspaces ===
  [fuel_centerline]  max a'dx = 2244.1   margin = -1372.4   too loose
  [t_avg_high_trip]  max a'dx =   33.3   margin =   -21.6   too loose
  [n_high_trip]      max a'dx = 1242.0   margin = -1242.0   too loose   (!!)
  [cold_leg_subcooled] max a'dx = 150.7  margin =  -135.7   too loose
\end{lstlisting}

The $n_{\mathrm{high\_trip}}$ row was the stopper: the barrier says
$n$ could deviate by $\pm 1242\times$ nominal.  Physically meaningless.
Fair enough to ask: is this because \emph{LQR isn't in the
barrier computation?}  The code uses $A_{\mathrm{cl}} = A - BK$, but
that detail was buried and Dane rightly flagged it.

\begin{decision}
\textbf{Write \texttt{barrier\_compare\_OL\_CL.m}} to explicitly
measure LQR's contribution to the barrier.  Same Qbar, same
disturbance envelope, two runs: open-loop $A$ vs.\ closed-loop
$A_{\mathrm{cl}} = A - BK$.
\end{decision}

Result (\cref{fig:ol-vs-cl}):

\begin{figure}[h]
  \centering
  \begin{tabular}{lrrr}
  \hline
  halfspace & open-loop & LQR closed-loop & ratio CL/OL \\
  \hline
  \texttt{fuel\_centerline}      & 26{,}941{,}862 K & 1{,}137 K & $4.2 \times 10^{-5}$ \\
  \texttt{t\_avg\_high\_trip}    & 788{,}220 K    & 33.2 K    & $4.2 \times 10^{-5}$ \\
  \texttt{n\_high\_trip}         & 27{,}379{,}931 & 1{,}242   & $4.5 \times 10^{-5}$ \\
  \texttt{cold\_leg\_subcooled}  & 1{,}806{,}166 K& 77.8 K    & $4.3 \times 10^{-5}$ \\
  $\gamma$ (ellipsoid level)     & $1.04 \times 10^{13}$ & $1.85 \times 10^{4}$ & $1.78 \times 10^{-9}$ \\
  \hline
  \end{tabular}
  \caption{Open-loop vs.\ LQR closed-loop Lyapunov barrier bounds on
  each \texttt{inv2\_holds} halfspace, same $\bar Q$ and same
  disturbance envelope.  LQR improves every bound by a uniform
  $\sim$20{,}000$\times$, but the bounds remain physically absurd.
  The open-loop version says $n$ could deviate by 27~million; the
  closed-loop says 1{,}242.  Both unusable.  The ceiling is plant
  anisotropy ($\Lambda/\tau_{\mathrm{thermal}} \sim 10^{-4}$ timescale
  ratio), not controller tuning.  Generated by
  \texttt{barrier\_compare\_OL\_CL.m}.}
  \label{fig:ol-vs-cl}
\end{figure}

\begin{derivation}
Why the eigenvalues barely move but $\gamma$ changes 9 orders of
magnitude:
\begin{itemize}
  \item $\max \Re(\mathrm{eig}\,A) = -0.0125$ — slowest thermal mode.
  \item $\max \Re(\mathrm{eig}\,A_{\mathrm{cl}}) = -0.0124$ — virtually
        identical.  LQR cannot speed up the slowest mode past what
        physics allows.
  \item But $\gamma = c_{\mathrm{inv}} \propto (w_{\mathrm{bar}} \sqrt{B_w^\top P B_w} / \mu)^2$,
        and $P$ changes \emph{shape} between OL and CL: LQR makes the
        disturbance penetrate $V$ far less.  $B_w^\top P B_w$ drops
        substantially under LQR, $\mu$ stays similar.  Ratio squared
        $\to$ nine orders of magnitude.
\end{itemize}
\end{derivation}

\begin{limitation}
Quadratic Lyapunov barrier is insufficient for this plant regardless
of LQR tuning.  Root cause: the plant's natural timescale spread
(prompt-neutron $\Lambda = 10^{-4}$~\unit{\second} vs.\ thermal $\sim 10$~\unit{\second})
forces any Lyapunov $P$ to be ill-conditioned by $\sim 10^4$.  The
resulting ellipsoid stretches absurdly in the fast directions when
projected to physical coordinates.  Fix:\ non-quadratic barrier
(polytopic or SOS).  Neither implemented yet.
\end{limitation}

\subsection*{Error 3: the heatup-rate invariant IS expressible as a halfspace}

In the first version of \texttt{predicates.json}, I put a comment on
\texttt{inv1\_holds} claiming a ramp-rate constraint could not be
written as a state halfspace without state augmentation.  That was
wrong.

\begin{derivation}
From \texttt{pke\_th\_rhs.m}, $\dot T_c$ is linear in
$(T_f, T_c, T_{\mathrm{cold}})$:
$$\dot T_c = \frac{hA(T_f - T_c) - 2 W c_c (T_c - T_{\mathrm{cold}})}{M_c c_c}
= a_f T_f + a_c T_c + a_{\mathrm{cold}} T_{\mathrm{cold}}$$
with
\begin{align*}
a_f     &= \frac{hA}{M_c c_c} = \frac{5 \times 10^7}{1.09 \times 10^8} = +0.4587~\unit{\per\second} \\
a_c     &= -\frac{hA + 2 W c_c}{M_c c_c} = -\frac{1.045 \times 10^8}{1.09 \times 10^8} = -0.9587~\unit{\per\second} \\
a_{\mathrm{cold}} &= \frac{2 W c_c}{M_c c_c} = \frac{5.45 \times 10^7}{1.09 \times 10^8} = +0.5000~\unit{\per\second}
\end{align*}
Coefficients sum to zero exactly: at thermal equilibrium
($T_f = T_c = T_{\mathrm{cold}}$), $\dot T_c = 0$.

A ramp-rate constraint $|\dot T_c| \leq r_{\max}$ is therefore two
halfspaces over the 10-state vector, with three nonzero coefficients
each (on indices $8, 9, 10$):
\begin{align*}
+a_f T_f + a_c T_c + a_{\mathrm{cold}} T_{\mathrm{cold}} &\leq r_{\max} \\
-a_f T_f - a_c T_c - a_{\mathrm{cold}} T_{\mathrm{cold}} &\leq r_{\max}
\end{align*}
Clean polyhedron; no augmentation.  The confusion came from conflating
this with rate constraints on \emph{non-linearly-driven} states — e.g.\
$|\dot n|$ involves $(\rho - \beta) n / \Lambda$, which is nonlinear,
and \emph{that} would need augmentation.  For the coolant temperature
rate, the thermal-hydraulic subsystem is linear-in-state and the
halfspace drops out directly.
\end{derivation}

$r_{\max}$ chosen as 50~\unit{\celsius\per\hour} (tech-spec nominal is
28~\unit{\celsius\per\hour}; add budget for transient overshoot).
\texttt{plant-model/test\_heatup\_rate.m} measures the actual peak
rate from the current \texttt{ctrl\_heatup.m} controller:

\begin{lstlisting}[style=terminal]
=== Heatup rate statistics ===
  max dT_c/dt:  0.0135 C/s  =   48.5 C/hr
  min dT_c/dt:  0.0000 C/s  =    0.0 C/hr
  rate limit:   +/- 0.01389 C/s = +/- 50.0 C/hr
  violates upper? false
  violates lower? false
\end{lstlisting}

Right at 50~\unit{\celsius\per\hour} peak during mid-ramp — barely
inside the budget, and \emph{70\% above tech-spec nominal}.  Would
fail a strict 28~\unit{\celsius\per\hour} interpretation.  Documented
as an open controller-tuning item.

\subsection*{The mode-obligation taxonomy}

The key conceptual move of the session.  Walking through the reach
code, Dane observed that the per-mode reach obligations aren't all the
same shape.  They split cleanly in two:

\begin{decision}
\textbf{Two kinds of mode, two kinds of obligation.}

\textbf{Equilibrium modes} (\texttt{q\_operation}, \texttt{q\_shutdown}):
plant is parked at a setpoint.  Obligation is forever-invariance under
bounded disturbance.  External disturbance matters — it's the thing
that could push the state out.  \emph{This is what the 2026-04-17
operation-mode reach actually proves} (up to the linearization
caveat).

\textbf{Transition modes} (\texttt{q\_heatup}, \texttt{q\_scram}):
plant is being driven from one mode's territory to another's.
Obligation is \emph{reach-avoid}: from $X_{\mathrm{entry}}$, reach
$X_{\mathrm{exit}}$ within bounded time $[T_{\min}, T_{\max}]$, while
maintaining $\mathrm{inv}_{\mathrm{mode}}$ along the way.  External
disturbance in these modes is essentially zero
($Q_{\mathrm{sg}} = 0$ during heatup, rods slam on scram).  The only
``disturbance'' is parametric uncertainty in the plant coefficients,
which we defer.
\end{decision}

The point of per-mode reach is \emph{not} generic disturbance
rejection — it's to prove that the DRC's discrete transitions
physically fire in finite time on the real plant.  The reach tube is
the artifact that transfers discrete correctness
(\texttt{ltlsynt}-guaranteed) to physical correctness.

Formal claim per transition mode (heatup example):
\begin{align}
&\text{Given } x(0) \in X_{\mathrm{entry}}(\text{heatup}), \notag \\
&\text{applying } u = \texttt{ctrl\_heatup}(t, x, \mathrm{plant}, \mathrm{ref}), \notag \\
&\text{with } Q_{\mathrm{sg}} = 0, \notag \\
&\text{we show } \exists\, T \in [T_{\min}, T_{\max}] \text{ s.t.} \notag \\
&\quad x(T) \in X_{\mathrm{exit}}(\text{heatup}) \quad \text{[reach]} \notag \\
&\quad x(t) \in \mathrm{inv1\_holds} \quad \forall t \in [0, T] \quad \text{[avoid]}
\end{align}

The entry/exit sets and time budgets need machine-readable
definitions.  Added a \texttt{mode\_boundaries} block to
\texttt{predicates.json}:
\begin{itemize}
  \item Each mode's \texttt{X\_entry\_polytope} as box bounds per state.
  \item Each mode's \texttt{X\_exit\_predicate} as a reference to named
        predicates.
  \item Each mode's \texttt{T\_min}, \texttt{T\_max} (\texttt{null} for
        equilibrium modes).
\end{itemize}

Values are engineering-reasonable placeholders:
\texttt{T\_max(heatup)} = 5~\unit{\hour} (tech-spec 28~\unit{\celsius/\hour}
over 60~\unit{\degreeFahrenheit} span is 1.19~\unit{\hour} nominal, 4$\times$
margin).  \texttt{T\_min(heatup)} = 2~\unit{\hour} 9~\unit{\minute} (any
faster violates the rate invariant).  \texttt{T\_max(scram)} = 60~\unit{\second}
(NRC expects seconds; 60 is generous for idealized rod insertion).

\begin{limitation}
All mode-boundary values are engineering guesses.  Thesis defense
will require real tech-spec numbers pinned to a specific plant.
Flagged in the JSON as \texttt{\_placeholder\_warning}.
\end{limitation}

\subsection*{WALKTHROUGH.md — standalone reach documentation}

With the mode-obligation taxonomy, predicate restructure, and barrier
findings in hand, wrote \texttt{reachability/WALKTHROUGH.md} as a
standalone document.  Captures the taxonomy table verbatim, per-mode
obligation definitions, code walkthrough, per-halfspace results, and
every known limitation.

\apass{WALKTHROUGH.md should probably be folded into a journal-style
chapter at some point, but for now it's the external-facing doc people
read first.  Keep it in sync when structural things change.}

\subsection*{Julia nonlinear reach — first attempt, partial success}

With linear work consolidated, turned to the real soundness question.
The linear reach proves the LQR closed-loop is safe, \emph{if} we
trust the linearization.  To close the gap, we need nonlinear reach.

Framework choice: \texttt{ReachabilityAnalysis.jl} with TMJets (Taylor
model integration).  Already have Julia scaffolding in
\texttt{julia-port/} and know the package.

Target: heatup mode reach, saturation disabled (simplifying
assumption to establish the framework; will add back saturation as
a hybrid sub-mode later).

\subsubsection*{The first three attempts failed}

\begin{deadend}
\textbf{Attempt 1:} naive RHS with \texttt{plant} as a function
argument.  Fails immediately with
\texttt{MethodError: setindex!(::Taylor1\{Float64\}, ::TaylorModel1\{...\}, ...)} —
Taylor model arithmetic needs the RHS in a specific form.  Need
\texttt{@taylorize}.

\textbf{Attempt 2:} add \texttt{@taylorize} to the function.  Same error.
Missing piece: the macro only handles \emph{constants} in arithmetic,
not struct-field accesses.  Need all plant parameters as \texttt{const}
globals.

\textbf{Attempt 3:} inline all constants.  Still fails.  Running out
of ideas; then noticed that \texttt{min()} inside the body (for the
ramp-reference clamp) is non-smooth — Taylor models can't handle
non-smooth operations.  Also: the raw time argument \texttt{t} in the
signature was interacting badly with TMJets' internal time parameter.
\end{deadend}

\begin{decision}
\textbf{Attempt 4 (landed):} three changes together.
\begin{enumerate}
  \item All plant parameters as \texttt{const} at module scope.
  \item \texttt{@taylorize} on the RHS, all arithmetic scalar (loops
        unrolled, no function calls in the hot path).
  \item Time carried as an \emph{augmented state} $x_{11}$ with
        $\dot x_{11} = 1$.  Instead of $T_{\mathrm{ref}}(t)$, the RHS
        references $T_{\mathrm{ref}}(x_{11})$ with no \texttt{min()}
        — valid while the ramp hasn't hit the target clamp, which is
        true for our probe horizons.
\end{enumerate}
\end{decision}

The final RHS body (see
\texttt{julia-port/scripts/reach\_heatup\_nonlinear.jl}):

\begin{lstlisting}[language=Julia,caption={Taylor-model-aware heatup RHS}]
@taylorize function rhs_heatup_taylor!(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] = one(x[1])     # d(time)/dt = 1, in the right Taylor type
    return nothing
end
\end{lstlisting}

Probes at $T = 10, 60, 300$~\unit{\second}:

\begin{lstlisting}[style=terminal]
--- Probe T = 10.0 s ---
  TMJets: 10583 reach-sets
  t reached:    10.0 s of 10.0
  n envelope:   [-0.0005178, 0.005015]
  T_f envelope: [274.46, 295.01] C
  T_c envelope: [274.45, 295.00] C
  T_cold env:   [270.00, 287.16] C

--- Probe T = 60.0 s ---
  Warning: Maximum number of integration steps reached; exiting.
  TMJets: 50000 reach-sets
  FAILED: AssertionError: radius must be nonnegative but is [NaN, NaN, ...]

--- Probe T = 300.0 s ---
  (same as 60s)
\end{lstlisting}

\textbf{10 seconds works}; 60 seconds onward exhaust the step budget
and then propagate NaN.  The 10-second envelope is sound — the
$n$-envelope going slightly negative is over-approximation tax
(physically impossible, numerically correct for a box hull).

\begin{derivation}
Step-size analysis.  TMJets' adaptive stepper shrinks $\Delta t$ to
resolve the fastest active mode.  Prompt-neutron timescale is
$\Lambda = 10^{-4}$~\unit{\second}.  For the Taylor remainder
(\texttt{abstol=1e-10}) to be bounded, the stepper needs
$\Delta t \lesssim 10^{-3}$~\unit{\second}.  Over a 10~\unit{\second}
horizon that's $\sim 10{,}000$ steps — consistent with the observed
10{,}583.

Extrapolate: a 5-hour heatup reach would need $\sim 1.8 \times 10^7$
steps at this cadence.  Totally infeasible, and we'd hit numerical
precision floor long before that.
\end{derivation}

\begin{limitation}
\textbf{TMJets on the full 10-state PKE is stuck at $\sim$10-second
horizons.}  Prompt-neutron stiffness ($\Lambda = 10^{-4}$~\unit{\second}
vs.\ thermal $\sim 10$~\unit{\second}, ratio $10^5$) forces tiny
time steps.  Step budget runs out, then numerical underflow
produces NaN in the precursor-decay arithmetic.  Not a framework
bug; not a tuning issue; a fundamental property of reach on the
full PKE.
\end{limitation}

\subsection*{The singular-perturbation remedy (not done yet)}

\begin{decision}
For hours-long reach, the prompt neutron timescale must be eliminated
from the state dynamics.  Standard move in reactor-kinetics reach
literature: set $\dot n \approx 0$ and solve for $n$ algebraically.
This gives the \emph{prompt-jump approximation} (sometimes written
as the \emph{reduced-order PKE}):
$$(\rho - \beta) n / \Lambda + \sum_i \lambda_i C_i = 0
\quad \Longrightarrow \quad
n \approx \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho}$$
valid when $\rho < \beta$ (subcritical).  $n$ becomes algebraic in
the remaining state; the dynamic state drops to 9-D, the stiffness
disappears, and TMJets can take steps on the thermal timescale
($\sim 1$~\unit{\second} per step is fine).

Tradeoff: we lose the prompt-supercritical regime, which is a
$\sim 50$~\unit{\microsecond} transient.  For heatup over hours, that
transient is irrelevant to safety.  Soundness is reduced (we no
longer capture prompt-dynamic trajectories), but the reduction
introduces a bounded error that can be characterized separately.
Deferred to a future session.
\end{decision}

\subsection*{Session result}

Artifacts:
\begin{itemize}
  \item \texttt{reachability/predicates.json} restructured:
        three-layer (deadbands, safety limits, mode invariants)
        + \texttt{mode\_boundaries} with per-mode entry/exit/time.
  \item \texttt{reachability/load\_predicates.m} handles the new format
        including multi-coefficient rows (for the rate halfspace).
  \item \texttt{reachability/reach\_operation.m} and
        \texttt{barrier\_lyapunov.m} now report per-halfspace margins
        against \texttt{inv2\_holds}.
  \item \texttt{reachability/barrier\_compare\_OL\_CL.m} — OL vs.\
        CL Lyapunov barrier.
  \item \texttt{reachability/WALKTHROUGH.md} — 550-line standalone
        document.
  \item \texttt{julia-port/scripts/reach\_heatup\_nonlinear.jl} and
        \texttt{sim\_heatup.jl}.  Nonlinear reach framework proven
        to 10~\unit{\second}.
  \item Six commits on \texttt{main}.
\end{itemize}

Key findings:
\begin{itemize}
  \item Quadratic Lyapunov barrier fails structurally on this plant.
        LQR improves it 20{,}000$\times$ but not enough.  Need
        polytopic or SOS barriers.
  \item Heatup rate is a clean state halfspace (three-coefficient row).
        Current controller peaks at 48.5~\unit{\celsius\per\hour} —
        tight against a strict 28-spec interpretation.
  \item Per-mode reach obligations split cleanly into equilibrium
        (forever-invariance under disturbance) vs.\ transition
        (reach-avoid with time budget, essentially zero disturbance).
  \item Nonlinear reach framework functional in Julia at 10-second
        horizons.  Prompt-neutron stiffness is the wall; reduced-order
        PKE is the known remedy.
\end{itemize}

\subsection*{Open at close}

\begin{itemize}
  \item Reduced-order PKE for long-horizon nonlinear reach.
  \item Heatup reach to $T_{\max} = 5$~\unit{\hour}.
  \item Scram reach (same reduction will apply).
  \item Add polytopic or SOS barrier as alternative certificate.
  \item Port MATLAB to Julia end-to-end; remove MATLAB dependency.
  \item Build a FRET-adjacent UI for exploring the predicates
        $\to$ halfspaces correspondence.
  \item Lab journal to document all of the above (this is what got
        done in the 2026-04-20 evening session — see next entry).
\end{itemize}