From 72143bcff55d92f2a0d6e5f10ea46d24a55f7774 Mon Sep 17 00:00:00 2001 From: Dane Sabo Date: Mon, 20 Apr 2026 16:16:53 -0400 Subject: [PATCH] predicates: add heatup rate invariant as a linear halfspace MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Earlier placeholder claimed ramp-rate limits weren't expressible as state halfspaces without augmentation. That was wrong: dT_c/dt is linear in (T_f, T_c, T_cold) directly from pke_th_rhs (no neutronics coupling), so |dT_c/dt| <= r_max is two clean halfspaces over x. Coefficients from pke_params: a_f = hA / (M_c*c_c) = +0.4587 /s a_c = -(hA + 2*W*c_c)/(M_c*c_c) = -0.9587 /s a_cold = 2*W*c_c / (M_c*c_c) = +0.5000 /s Sum = 0 exact (equilibrium when all T's equal). Limit chosen: +/- 50 C/hr (tech-spec 28 C/hr + transient overshoot budget). Verified on actual heatup sim: max dT_c/dt = 48.5 C/hr, min = 0 C/hr. Passes our placeholder but tight — a strict 28 C/hr tech- spec invariant would be violated by current ctrl_heatup tuning (overshoot factor ~1.7x during mid-ramp). Generalized load_predicates.m to accept multi-coefficient halfspace rows via "row": [[state_idx, coeff], ...] format, in addition to the existing single-coefficient {state_index, coeff} form. Backward compatible. inv1_holds now conjoins fuel_centerline, cold_leg_subcooled, and the two rate halfspaces. DNBR still not modeled (would need an augmented predicate with a correlation-based safety margin). Hacker-Split: Dane asked about heatup rate invariant; realizing my earlier 'needs state augmentation' claim was wrong and the rate constraint is already linear. Fix it, verify against actual sim. Co-Authored-By: Claude Opus 4.7 (1M context) --- plant-model/test_heatup_rate.m | 20 ++++++++++++++++++++ reachability/load_predicates.m | 13 ++++++++++++- reachability/predicates.json | 23 +++++++++++++++++++++-- 3 files changed, 53 insertions(+), 3 deletions(-) create mode 100644 plant-model/test_heatup_rate.m diff --git a/plant-model/test_heatup_rate.m b/plant-model/test_heatup_rate.m new file mode 100644 index 0000000..9978d72 --- /dev/null +++ b/plant-model/test_heatup_rate.m @@ -0,0 +1,20 @@ +%% test_heatup_rate.m — measure the actual heatup-rate trajectory from ctrl_heatup +clear; addpath('controllers'); addpath('../reachability'); +plant = pke_params(); +pred = load_predicates(plant); +x0 = [1e-3; (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * 1e-3; ... + pred.constants.T_standby; pred.constants.T_standby; pred.constants.T_standby]; +ref = struct('T_start', pred.constants.T_standby, ... + 'T_target', plant.T_c0, ... + 'ramp_rate', 28/3600); +[t, X, U] = pke_solver(plant, @(t) 0, @ctrl_heatup, ref, [0 5400], x0); + +a_f = 0.4587; a_c = -0.9587; a_cold = 0.5000; +dTcdt = a_f*X(:,8) + a_c*X(:,9) + a_cold*X(:,10); % C/s + +fprintf('\n=== Heatup rate statistics ===\n'); +fprintf(' max dT_c/dt: %.4f C/s = %6.1f C/hr\n', max(dTcdt), max(dTcdt)*3600); +fprintf(' min dT_c/dt: %.4f C/s = %6.1f C/hr\n', min(dTcdt), min(dTcdt)*3600); +fprintf(' rate limit: +/- 0.01389 C/s = +/- 50.0 C/hr\n'); +fprintf(' violates upper? %s\n', string(any(dTcdt > 0.01389))); +fprintf(' violates lower? %s\n', string(any(dTcdt < -0.01389))); diff --git a/reachability/load_predicates.m b/reachability/load_predicates.m index 1611c0d..f5c8e5c 100644 --- a/reachability/load_predicates.m +++ b/reachability/load_predicates.m @@ -93,7 +93,18 @@ function group_out = parse_group(group_in, T_c0, T_f0, T_cold0, T_standby) b_poly = zeros(n_hs, 1); for i = 1:n_hs hs = get_hs(i); - A_poly(i, hs.state_index) = hs.coeff; + if isfield(hs, 'row') + % Multi-coefficient halfspace: hs.row is a k x 2 matrix + % where each row is [state_index, coeff] (jsondecode from + % JSON array-of-arrays). + coeffs = hs.row; + for r = 1:size(coeffs, 1) + A_poly(i, coeffs(r, 1)) = coeffs(r, 2); + end + else + % Single-coefficient halfspace: hs.state_index + hs.coeff. + A_poly(i, hs.state_index) = hs.coeff; + end b_poly(i) = evalin_context(hs.rhs_expr, T_c0, T_f0, T_cold0, T_standby); end group_out.(name).A_poly = A_poly; diff --git a/reachability/predicates.json b/reachability/predicates.json index 30e46b8..ea569ed 100644 --- a/reachability/predicates.json +++ b/reachability/predicates.json @@ -110,6 +110,20 @@ "halfspaces": [ { "state_index": 10, "coeff": 1.0, "rhs_expr": "T_cold0 + 15.0" } ] + }, + "heatup_rate_upper": { + "meaning": "Coolant heatup rate does not exceed tech-spec limit + overshoot margin.", + "concretization": "dT_c/dt = a_f*T_f + a_c*T_c + a_cold*T_cold <= 0.01389 C/s (50 C/hr; tech-spec 28 C/hr + transient overshoot budget)", + "_derivation": "dT_c/dt is linear in (T_f, T_c, T_cold) from pke_th_rhs.m: a_f=hA/(M_c*c_c)=+0.4587/s, a_c=-(hA+2*W*c_c)/(M_c*c_c)=-0.9587/s, a_cold=2*W*c_c/(M_c*c_c)=+0.5000/s. Coefficients sum to zero by construction (equilibrium when all T's equal).", + "halfspaces": [ + { "row": [[8, 0.4587], [9, -0.9587], [10, 0.5000]], "rhs_expr": "0.01389" } + ] + }, + "heatup_rate_lower": { + "meaning": "Coolant cooldown rate during heatup doesn't exceed -50 C/hr (flag runaway cooldown).", + "halfspaces": [ + { "row": [[8, -0.4587], [9, 0.9587], [10, -0.5000]], "rhs_expr": "0.01389" } + ] } }, @@ -120,8 +134,13 @@ ], "inv1_holds": { "meaning": "Heatup mode safety envelope.", - "conjunction_of": ["fuel_centerline", "cold_leg_subcooled"], - "_placeholder_note": "Ramp-rate limit and DNBR not expressible as state halfspaces without augmentation." + "conjunction_of": [ + "fuel_centerline", + "cold_leg_subcooled", + "heatup_rate_upper", + "heatup_rate_lower" + ], + "_note": "dT_c/dt is linear in (T_f, T_c, T_cold), so ramp-rate halfspace has 3 nonzero coefficients per row. DNBR not modeled — would need an augmented correlation-based predicate." }, "inv2_holds": { "meaning": "Operation mode safety envelope.",