8d2c7d0956
Morning-review point 3 result: tight-entry heatup PJ reach with
Q_sg in [0, 5% P0] as a bounded parameter (augmented state x[10]).
T=60s: 7042 sets in 394s — T_c in [270.97, 291.0] — low-trip × loose
T=300s: 100k sets budget exhausted in 5400s —
T_c in [219.4, 316.3] — low-trip × loose
Compared to the no-disturbance tight-entry run (all 6 halfspaces at
300s, T_c in [281.05, 291.0]), the bounded steam-dump disturbance
costs the low-T_avg-trip discharge even at 60s. Physically correct
— steam dump pulls heat through secondary, cascades into cold-leg
and T_avg. The reach tube accurately captures this coupling.
Thesis-relevant finding: without modeled disturbance bounds, reach
tubes are over-optimistic. Quantifies how much of the prior
"all 6 halfspaces" result came from Q_sg=0 simplification vs.
actual controller robustness.
Results saved to results/reach_heatup_pj_with_steam_dump.mat.
Journal entry updated with the per-horizon table + decision box on
what this means for thesis claims.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
268 lines
13 KiB
TeX
268 lines
13 KiB
TeX
% ---------------------------------------------------------------------------
|
||
% 2026-04-21 — Polytopic & SOS barriers; Tikhonov bound for prompt-jump
|
||
% Live / B-style entry, A-style on Tikhonov derivation.
|
||
% ---------------------------------------------------------------------------
|
||
|
||
\session{2026-04-21 (overnight cont.)}{autonomous while Dane is at the gym}{%
|
||
Explore polytopic and SOS polynomial barriers on the operation mode,
|
||
work out the Tikhonov singular-perturbation bound that would make the
|
||
prompt-jump reduction rigorous rather than empirical, and leave
|
||
everything committed and documented.}
|
||
|
||
\section{2026-04-21 --- Polytopic / SOS barriers + Tikhonov bound}
|
||
\label{sec:20260421-polytopic-sos-tikhonov}
|
||
|
||
\subsection*{Polytopic barrier: naive check, expected failure mode}
|
||
|
||
Wrote \texttt{code/scripts/barrier\_polytopic.jl}. Test: is the polytope
|
||
$P = \texttt{inv2\_holds} \cap (\text{precursor tube-bounds})$ forward-invariant
|
||
under $A_{\mathrm{cl}}$ closed-loop with bounded $Q_{\mathrm{sg}}$?
|
||
Nagumo's theorem: check each face $a_i^\top x = b_i$ with an LP asking
|
||
for $\max\ a_i^\top (A_{\mathrm{cl}} x + B_w w)$ over the polytope and
|
||
admissible $w$. If $\leq 0$ for every face, $P$ is invariant.
|
||
|
||
Result: $2 / 18$ faces pass. The other $16$ can be crossed: LQR can't
|
||
contract a point on the safety boundary back inward because the
|
||
polytope includes regions far outside what the LQR reach can actually
|
||
reach. \textbf{Expected:} safety halfspaces + reach-tube-bounds
|
||
together form a set much larger than the actual minimal invariant,
|
||
so local outward velocities are plentiful.
|
||
|
||
\begin{decision}
|
||
The right approach for a tight polytopic barrier is \textbf{Blanchini's
|
||
pre-image algorithm}: $P_{k+1} = P_k \cap \{x : A_{\mathrm{cl}} x + B_w w \in P_k\ \forall w \in W\}$,
|
||
iterating until fixed point. The fixed point is the \emph{maximal
|
||
robustly controllable invariant set} inside $P_0 = $ safety polytope.
|
||
Each iteration adds faces; polytope combinatorial complexity grows.
|
||
Requires \texttt{Polyhedra.jl} + \texttt{CDDLib} for polytope ops,
|
||
HiGHS for LPs. 2--3 days focused work. Deferred.
|
||
\end{decision}
|
||
|
||
The naive check is not a failure; it's a diagnostic that tells us which
|
||
algorithmic tool we actually need.
|
||
|
||
\subsection*{SOS polynomial barrier: first success}
|
||
|
||
Wrote \texttt{code/scripts/barrier\_sos\_2d.jl}. Use \texttt{SumOfSquares.jl}
|
||
+ CSDP to find a polynomial $B(x)$ satisfying the Prajna--Jadbabaie
|
||
conditions:
|
||
\begin{enumerate}
|
||
\item $B(x) \leq 0$ on $X_{\mathrm{entry}}$.
|
||
\item $B(x) \geq 0$ on $X_{\mathrm{unsafe}}$ (complement of safety).
|
||
\item $\nabla B \cdot f \leq 0$ on $\{B = 0\}$.
|
||
\end{enumerate}
|
||
|
||
Reduced the operation-mode problem to a 2-state projection $(\delta n,
|
||
\delta T_c)$ after LQR, dropping the other 8 states (and therefore the
|
||
disturbance coupling, since $B_w$ projects to zero in this subset). Set
|
||
safety $|\delta T_c| \leq 5$~\unit{\kelvin} and $|\delta n| \leq 0.15$,
|
||
entry $|\delta T_c| \leq 0.1$ and $|\delta n| \leq 0.01$, unsafe
|
||
$\delta n \geq 0.15$ (high-flux-trip direction).
|
||
|
||
Technical simplification: instead of the bilinear Putinar form
|
||
$-(\nabla B \cdot f) - \sigma_b \cdot B$ SOS (which requires iterative
|
||
BMI decomposition), used the stronger condition $-(\nabla B \cdot f)$
|
||
SOS globally. Safe for linear Hurwitz closed-loop because such
|
||
systems admit a decreasing Lyapunov-like polynomial everywhere.
|
||
|
||
\textbf{Result:} CSDP returned \texttt{OPTIMAL}. A degree-4 polynomial
|
||
barrier exists:
|
||
|
||
\begin{lstlisting}[style=terminal, breaklines=true]
|
||
B(x) = -0.7596 + 15.149*x2^2 + 0.5816*x1*x2 + 35.2614*x1^2
|
||
- 0.1618*x2^3 + 7.0328*x1*x2^2 - 0.1035*x1^2*x2
|
||
+ 15.8024*x1^3 + 46.8212*x2^4 - 0.0107*x1*x2^3
|
||
+ 6.5748*x1^2*x2^2 - 0.1111*x1^3*x2 + 5.9248*x1^4
|
||
\end{lstlisting}
|
||
|
||
where $x_1 = \delta n$, $x_2 = \delta T_c$. Constant term negative
|
||
(\emph{B} at origin is negative, origin is in entry set); quartic in
|
||
$x_1$ dominates when $|\delta n|$ is large (pushing $B$ positive at
|
||
unsafe). \textbf{First non-quadratic barrier certificate for this
|
||
plant.}
|
||
|
||
\begin{limitation}
|
||
2D projection loses the precursor--thermal coupling and the disturbance
|
||
(which only enters $T_{\mathrm{cold}}$, projected out). Not a direct
|
||
safety claim for the 10-state system. Scaling to the full 10 states:
|
||
degree-4 monomials in 10 variables is $\binom{14}{4} = 1001$; the SDP
|
||
matrix is $\sim 1000 \times 1000$, which CSDP may struggle with.
|
||
Switching to Mosek (if licensed) or SCS (open source) would help.
|
||
The Putinar boundary form is the right long-term formulation;
|
||
iterative BMI solvers (PENBMI, iterative SOS) are the path.
|
||
\end{limitation}
|
||
|
||
\apass{Extend to full 10-state, keep degree 4 or reduce to degree 3,
|
||
add disturbance (via Schur complement or worst-case polytopic
|
||
bound), and iterate the Putinar/BMI solver until convergence. Probably
|
||
a week of focused work once the approach is chosen.}
|
||
|
||
\subsection*{Tikhonov bound for the prompt-jump reduction}
|
||
|
||
\begin{derivation}
|
||
Write the 10-state PKE in standard singular-perturbation form. Let
|
||
$y = n$ (fast) and $x = [C_1, \ldots, C_6, T_f, T_c, T_{\mathrm{cold}}]^\top$
|
||
(slow). The neutron balance is
|
||
$$\dot y = \frac{\rho(x) - \beta}{\Lambda}\, y + \sum_i \lambda_i C_i.$$
|
||
Multiplying through by $\Lambda$:
|
||
$$\Lambda \dot y = -(\beta - \rho(x)) y + \Lambda \sum_i \lambda_i C_i.$$
|
||
With $\varepsilon := \Lambda$ as the small parameter, and defining
|
||
$$g(x, y) := -(\beta - \rho(x))\, y + \varepsilon \sum_i \lambda_i C_i,$$
|
||
the system is
|
||
$$\dot x = f(x, y), \qquad \varepsilon \dot y = g(x, y),$$
|
||
exactly the form for \textbf{Tikhonov's theorem}.
|
||
|
||
The quasi-steady-state manifold is $g(x, y) = 0$:
|
||
$$y = h(x) := \frac{\varepsilon \sum_i \lambda_i C_i}{\beta - \rho(x)}
|
||
= \frac{\Lambda \sum_i \lambda_i C_i}{\beta - \rho(x)}.$$
|
||
This is exactly our prompt-jump formula for $n_{\mathrm{PJ}}$.
|
||
|
||
\textbf{Asymptotic stability of the fast subsystem} (with $x$ frozen):
|
||
$\frac{d(y - h(x))}{d\tau} = -(\beta - \rho(x))(y - h(x)) / \varepsilon$,
|
||
using $\tau = t/\varepsilon$ (fast time). Decay rate $(\beta - \rho)/\varepsilon$.
|
||
Since $\beta - \rho > 0$ (by the \texttt{prompt\_critical\_margin\_heatup}
|
||
invariant, conjoined into \texttt{inv1\_holds} as of this morning), the
|
||
fast dynamics are exponentially stable with time constant
|
||
$\varepsilon / (\beta - \rho) \leq \Lambda / (0.5\beta) \approx 3 \times 10^{-2}~\unit{\second}$.
|
||
|
||
\textbf{Tikhonov's theorem} (Khalil, \emph{Nonlinear Systems}, Thm 11.1;
|
||
Kokotović, Khalil, \& O'Reilly \emph{Singular Perturbation Methods in
|
||
Control}): under the hypotheses above, for sufficiently small $\varepsilon > 0$
|
||
and on any compact time interval $[t_1, T]$ after the boundary layer,
|
||
there exist positive constants $K_1, K_2$ such that
|
||
\begin{align}
|
||
|y(t) - h(\bar x(t))| &\leq K_1 \cdot \varepsilon + K_2 \cdot e^{-\gamma t / \varepsilon}, \\
|
||
|x(t) - \bar x(t)| &\leq K_3 \cdot \varepsilon,
|
||
\end{align}
|
||
where $\bar x$ is the reduced-system solution and $\gamma$ is the
|
||
fast-subsystem decay rate. After the initial layer $O(\varepsilon \log(1/\varepsilon))$,
|
||
the second term decays below the first and the error is uniformly
|
||
$O(\varepsilon) = O(\Lambda)$.
|
||
|
||
\textbf{Sanity check against our empirical validation.} With
|
||
$\Lambda = 10^{-4}$~\unit{\second} and typical problem magnitudes:
|
||
\begin{itemize}
|
||
\item Absolute error on $n$: $|n(t) - n_{\mathrm{PJ}}(t)| \leq K_1 \cdot 10^{-4}$
|
||
for some constant $K_1$. Our empirical max at $t = 1200$~\unit{\second}
|
||
was $|3.414 \times 10^{-3} - 3.410 \times 10^{-3}| \approx 4 \times 10^{-6}$.
|
||
If $K_1 \approx 40$, the bound is $4 \times 10^{-3}$; our data sits
|
||
three orders of magnitude tighter, consistent with $K_1$ being
|
||
plant-dependent and the actual error being substantially below
|
||
the worst-case bound.
|
||
\item Absolute error on temperatures: $|T(t) - \bar T(t)| \leq K_3 \cdot 10^{-4}$.
|
||
Empirical max was $7 \times 10^{-3}$~\unit{\kelvin}. If $K_3 \approx 70$,
|
||
this is consistent.
|
||
\end{itemize}
|
||
The constants $K_1, K_3$ are problem-dependent and bounded on the
|
||
reach set. A tight numerical estimate would require computing the
|
||
Jacobians of $f$ and $h$ along the trajectory; rough back-of-envelope
|
||
from the empirical data gives the bound meaningful physical interpretation.
|
||
\end{derivation}
|
||
|
||
\begin{decision}
|
||
\textbf{For the thesis:} state the PJ error as
|
||
$\|x(t) - x_{\mathrm{PJ}}(t)\| \leq C \Lambda = O(10^{-4})$
|
||
\emph{in state units}, invoking Tikhonov's theorem with the
|
||
\texttt{prompt\_critical\_margin\_heatup} invariant (proven by
|
||
reach) as the hypothesis. The constant $C$ can be bounded above by
|
||
problem-specific norms of the Jacobians of $f, h$ restricted to the
|
||
reach set, which are themselves polytope-bounded state functions
|
||
and thus computable.
|
||
|
||
This upgrades the validation-based ``we ran it and 0.1\% was the max''
|
||
to a rigorous ``bounded by $C \Lambda$ where $C$ depends on properties
|
||
of the reach set, themselves bounded by the safety halfspaces.''
|
||
|
||
\textbf{Remaining gap}: compute $C$ numerically on our reach tube.
|
||
Straightforward: evaluate $\partial f / \partial y$ and $\partial h / \partial x$
|
||
at the vertices of $X_{\mathrm{entry}}$ + reach envelope, take the max.
|
||
One-session task.
|
||
\end{decision}
|
||
|
||
\subsection*{Other odds and ends}
|
||
|
||
\textbf{Scram entry-set expansion (user's morning point 2) ---
|
||
landed later this session.} Built
|
||
\texttt{code/scripts/reach/reach\_loca\_operation.jl} (LQR reach under
|
||
$Q_{\mathrm{sg}} \in [0,\ 1.5 P_0]$, the steam-line-break envelope,
|
||
for 3~\unit{\second} horizon) and
|
||
\texttt{code/scripts/reach/reach\_scram\_pj\_fat.jl} (bounding-box
|
||
union of hot-standby + heatup-tight envelope + operation-LQR envelope
|
||
+ LOCA envelope, clamps LOCA's numerical outliers on precursors to
|
||
physically plausible bounds, reruns PJ scram reach).
|
||
|
||
Result: $n$ decays monotonically
|
||
\textbf{$0.047 \to 0.021 \to 0.0094$ over $\{10, 30, 60\}$~\unit{\second}},
|
||
factor-of-five per minute even starting from the fat entry (which
|
||
includes the LOCA-perturbed post-operation state). Temperatures
|
||
fall from $[226,\ 361]\,^\circ\mathrm{C}$ (clamp-saturated initially)
|
||
toward the decay-heat equilibrium. No step-budget truncation;
|
||
23{,}919 reach-sets over 60~\unit{\second}. $X_{\mathrm{exit}}$
|
||
threshold of $n \leq 10^{-4}$ still not reached in 60~\unit{\second};
|
||
same $T_{\max}$-vs-plant-decay mismatch flagged in the earlier scram
|
||
entry. \emph{But}: this is now a defensible obligation because the
|
||
entry set represents ``anywhere the plant could realistically be,''
|
||
not the narrow 1~K box around $x_{\mathrm{op}}$ of the earlier run.
|
||
|
||
\apass{The LOCA reach itself is numerically loose (box-hull propagation
|
||
amplifies slow precursor modes under large disturbance, so horizons
|
||
$> 3$~\unit{\second} blow up). A proper zonotope-generator
|
||
propagator would fix this; deferred.}
|
||
|
||
\textbf{Heatup with steam-dump $Q_{\mathrm{sg}}$ demand (user's morning
|
||
point 3) --- landed this session.} Built
|
||
\texttt{code/configs/heatup/with\_steam\_dump.toml} +
|
||
\texttt{code/scripts/reach/reach\_heatup\_pj\_sd.jl}: 11-state RHS
|
||
(9 physics + $x_{10} = Q_{\mathrm{sg}}$ as augmented bounded parameter
|
||
with $\dot x_{10} = 0$, $x_{11} = t$). Entry box on $Q_{\mathrm{sg}}$:
|
||
$[0,\ 0.05 P_0]$ (steam dump to atmosphere, conservative).
|
||
|
||
Results from the tight X\_entry + steam-dump run:
|
||
\begin{lstlisting}[style=terminal]
|
||
--- Probe T = 60.0 s ---
|
||
TMJets: 7042 reach-sets in 393.6 s
|
||
T_c envelope: [270.97, 291.0] °C
|
||
Low-T_avg trip (T_c >= 280): × loose
|
||
|
||
--- Probe T = 300.0 s ---
|
||
Max-steps budget exhausted (100,000 reach-sets, 5403 s wall)
|
||
T_c envelope: [219.4, 316.28] °C
|
||
Low-T_avg trip: × loose
|
||
\end{lstlisting}
|
||
|
||
\textbf{Steam-dump disturbance costs the low-$T_{\mathrm{avg}}$ trip
|
||
discharge even at 60~\unit{\second}.} Without the dump
|
||
($Q_{\mathrm{sg}} = 0$ exact), the tight-entry run cleared all six
|
||
halfspaces at 300~\unit{\second} with T\_c $\in [281.05, 291.0]$. With
|
||
the dump in $[0, 5\%]$, T\_c lower bound drops to 270.97~$^\circ$C ---
|
||
physically consistent: steam dump pulls heat from secondary, cools
|
||
cold-leg, cascades into T\_avg.
|
||
|
||
At 300~\unit{\second} with the dump, step budget exhausts (100k sets
|
||
in 90~\unit{\minute} wall) and the envelope blows out. Bigger budget
|
||
or entry-box refinement would likely recover; deferred.
|
||
|
||
\begin{decision}
|
||
The steam-dump result is pedagogically useful for the thesis: it
|
||
shows quantitatively how much of the safety margin comes from
|
||
``plant is isolated'' modeling vs.\ realistic operational
|
||
assumptions. Without accurate disturbance bounds the reach tube
|
||
is over-optimistic.
|
||
\end{decision}
|
||
|
||
\apass{The reach tube plots (Dane's point 4) for the heatup PJ tight
|
||
entry revealed a controller-reference mismatch: with
|
||
$X_{\mathrm{entry}}$ at $T_c \in [285, 291]$ and the controller's
|
||
ramp reference starting at $T_{\mathrm{standby}} = 275$, the
|
||
feedback-lin controller commands cooling ($\rho < 0$ throughout the tube).
|
||
The heatup physics isn't captured. Fix: parameterize the controller's
|
||
\texttt{T\_start} from the current $T_c$ at mode entry. Documented
|
||
in the tube-plot commit message.}
|
||
|
||
\subsection*{Remote push blocked, commits all local}
|
||
|
||
The harness correctly blocked an agent-inferred gitea URL when I tried
|
||
to push for backup. Flagged in \texttt{OVERNIGHT\_NOTES.md} with the
|
||
exact command Dane needs to run. All work is committed locally on
|
||
\texttt{main}; nothing lost.
|