From a20d2a05e91d3f2330425900a7a8ab924f1e0576 Mon Sep 17 00:00:00 2001 From: Dane Sabo Date: Mon, 20 Apr 2026 16:04:40 -0400 Subject: [PATCH] predicates: split operational deadbands from hard safety limits MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously conflated two different kinds of constraint: - operational deadbands (|T_c - T_c0| <= 5 F) used by the DRC for mode transitions. Symmetric bands around setpoint. Violating these is an operator/operational issue, not a safety issue. - safety limits (T_f <= 1200 C, T_c <= 320 C, n <= 1.15, etc.) are hard one-sided halfspaces corresponding to physical damage mechanisms or reactor-trip setpoints. THESE are what a safety barrier/reach must discharge. predicates.json now has three groups: - operational_deadbands (t_avg_above_min, t_avg_in_range, p_above_crit) - safety_limits (fuel_centerline, t_avg_high_trip, t_avg_low_trip, n_high_trip, n_low_operation, cold_leg_subcooled) - mode_invariants (inv1_holds, inv2_holds as conjunctions of safety_limits) reach_operation.m and barrier_lyapunov.m both now report halfspace-by- halfspace margins against inv2_holds. Attributable failure analysis: we can see WHICH limit is tightest. Reach tube (under +/-15% Q_sg load): passes all 6 safety halfspaces. Tightest margin is n_high_trip at +0.138 (12% from trip). Temperature directions have 10-870 K margin. Lyapunov barrier (same): fails all 6. Worst is n_high_trip with -2365 margin — the ellipsoid says n could deviate by +/-2364, which is physically meaningless. Anisotropy cost made visible per-direction. Motivates SOS / polytopic barriers for the thesis chapter. load_predicates.m now returns .operational_deadbands, .safety_limits, and .mode_invariants. Existing callers that only used .constants or .t_avg_in_range still work because those live under the old keys. Hacker-Split: user caught that the barrier was checking the wrong invariant; safety limits != operating deadband. Restructured so the proof target matches the physical claim. Co-Authored-By: Claude Opus 4.7 (1M context) --- docs/figures/reach_operation_Tc.png | Bin 63233 -> 63233 bytes docs/figures/reach_operation_n.png | Bin 35326 -> 35326 bytes reachability/barrier_lyapunov.m | 34 +++++- reachability/load_predicates.m | 94 +++++++++++----- reachability/predicates.json | 163 +++++++++++++++++++--------- reachability/reach_operation.m | 30 ++++- 6 files changed, 233 insertions(+), 88 deletions(-) diff --git a/docs/figures/reach_operation_Tc.png b/docs/figures/reach_operation_Tc.png index ee90c364039c5ba60a33b135ce4e387e84c53cbb..9fdc050bc80dc51f48c707dd3452ba2c13101488 100644 GIT binary patch delta 35 rcmZpC$K3did4dzW2s8ILh0@g%J%hQ;tPG5;3=O%%Z_nAdJn9<&<1Gz7 delta 35 rcmZpC$K3did4dzWFo)r_>B->}J%hPTtqd%zOf2UKS(R^G9`y|Xr&geJiHSC-p&j+ delta 35 rcmex2nd#qTrU_2$!W>5a^Nik3^bF=UwKA} 0), 'P not positive definite'); %% ===== Safety spec (used by sweep and final check) ===== +% Load inv2_holds (conjunction of safety halfspaces) from the predicates +% source of truth. Each row k of A_inv defines a halfspace a_k' x <= b_k; +% the barrier must bound max(a_k' * dx) over the ellipsoid for each k. +addpath('../plant-model'); +pred = load_predicates(plant); +inv2 = pred.mode_invariants.inv2_holds; +A_inv = inv2.A_poly; +b_inv = inv2.b_poly; +comp_names = inv2.components; +n_halfspaces = size(A_inv, 1); + +% Convert limits to deviation from x_op: +% for halfspace a' x <= b, the deviation-frame bound is a' dx <= b - a' x_op. +b_inv_dev = b_inv - A_inv * x_op; + +% Backward-compat scalars for existing prints. e9 = zeros(10, 1); e9(9) = 1; delta_safe_Tc = 5.0; @@ -161,8 +177,24 @@ else fprintf(' Try a tighter LQR (bigger Q_Tc or smaller R) or tighter X_entry.\n'); end +%% ===== Per-halfspace barrier check against inv2_holds ===== +% For each safety halfspace a' dx <= b_dev, the max of a' dx over the +% gamma-ellipsoid {dx : dx' P dx <= gamma} is sqrt(gamma * a' P^{-1} a). +% Compare to b_dev (the headroom from x_op to the safety limit). +fprintf('\n=== Lyapunov barrier vs inv2_holds halfspaces ===\n'); +Pinv = inv(P); +for k = 1:n_halfspaces + a = A_inv(k, :).'; + max_adx = sqrt(gamma * (a.' * Pinv * a)); + margin = b_inv_dev(k) - max_adx; + status = 'OK'; + if margin < 0, status = '*** BARRIER TOO LOOSE ***'; end + fprintf(' [%-20s] headroom = %8.3f | max a''dx = %8.3f | margin = %+8.3f %s\n', ... + comp_names{k}, b_inv_dev(k), max_adx, margin, status); +end + save(fullfile('.', 'barrier_lyapunov_result.mat'), ... 'P', 'gamma', 'c_entry', 'c_inv', 'w_bar', 'K', 'A_cl', 'delta_entry', ... - 'max_dTc_on_ellipsoid', 'delta_safe_Tc', '-v7'); + 'max_dTc_on_ellipsoid', 'delta_safe_Tc', 'A_inv', 'b_inv_dev', '-v7'); fprintf('\nSaved barrier to ./barrier_lyapunov_result.mat\n'); diff --git a/reachability/load_predicates.m b/reachability/load_predicates.m index 79bc6b0..1611c0d 100644 --- a/reachability/load_predicates.m +++ b/reachability/load_predicates.m @@ -1,25 +1,21 @@ function pred = load_predicates(plant) % LOAD_PREDICATES Read predicates.json and resolve rhs_expr into numbers. % -% Each halfspace entry in the JSON stores rhs_expr as a string because -% several of the bounds are defined relative to plant-derived constants -% (T_c0, T_cold0, T_standby). We evaluate those expressions here in a -% controlled workspace that exposes exactly those names plus the -% derived offsets from the JSON. +% Returns: +% pred.constants - struct with T_c0, T_cold0, T_f0, T_standby +% pred.operational_deadbands - struct of predicates (each: A_poly, b_poly, meaning) +% pred.safety_limits - struct of halfspace limits (each: A_poly, b_poly, meaning) +% pred.mode_invariants - struct mapping mode-invariant name to +% its conjoined (A_poly, b_poly, components) % -% Returns a struct `pred` with one field per predicate, each a struct -% holding an n_hs x 10 matrix A_poly and an n_hs x 1 vector b_poly such -% that the predicate is { x : A_poly * x <= b_poly }. -% -% Also returns: -% pred.constants - a struct with T_c0, T_cold0, T_f0, T_standby, etc. +% Each halfspace is of the form { x : A_poly * x <= b_poly }. Conjunctions +% (polytopes) stack rows into A_poly / b_poly. % % Usage: % plant = pke_params(); % pred = load_predicates(plant); -% A = pred.t_avg_in_range.A_poly; % 2 x 10 -% b = pred.t_avg_in_range.b_poly; % 2 x 1 -% is_in = all(A * x <= b); % predicate check on state x +% inv2 = pred.mode_invariants.inv2_holds; +% is_in = all(inv2.A_poly * x <= inv2.b_poly); here = fileparts(mfilename('fullpath')); raw = fileread(fullfile(here, 'predicates.json')); @@ -33,23 +29,59 @@ function pred = load_predicates(plant) T_standby = T_c0 + T_standby_offset_C; %#ok pred.constants = struct( ... - 'T_c0', plant.T_c0, ... - 'T_f0', plant.T_f0, ... - 'T_cold0', plant.T_cold0, ... - 'T_standby', T_standby, ... + 'T_c0', plant.T_c0, ... + 'T_f0', plant.T_f0, ... + 'T_cold0', plant.T_cold0, ... + 'T_standby', T_standby, ... 'T_standby_offset_C', T_standby_offset_C, ... 'T_standby_offset_F', J.derived.T_standby_offset_F, ... 't_avg_in_range_halfwidth_C', J.derived.t_avg_in_range_halfwidth_C, ... - 'p_above_crit_threshold_n', J.derived.p_above_crit_threshold_n); + 'p_above_crit_threshold_n', J.derived.p_above_crit_threshold_n, ... + 'T_fuel_limit_C', J.derived.T_fuel_limit_C, ... + 'T_c_high_trip_C', J.derived.T_c_high_trip_C, ... + 'n_high_trip', J.derived.n_high_trip); - % --- Loop over predicates, build A/b matrices --- - names = fieldnames(J.predicates); + % --- operational_deadbands --- + pred.operational_deadbands = parse_group(J.operational_deadbands, ... + T_c0, T_f0, T_cold0, T_standby); + + % --- safety_limits --- + pred.safety_limits = parse_group(J.safety_limits, ... + T_c0, T_f0, T_cold0, T_standby); + + % --- mode_invariants: conjunctions of safety_limits entries --- + inv_names = fieldnames(J.mode_invariants); + for k = 1:numel(inv_names) + name = inv_names{k}; + if startsWith(name, '_') || startsWith(name, 'x_'), continue, end + entry = J.mode_invariants.(name); + if ~isstruct(entry) || ~isfield(entry, 'conjunction_of'), continue, end + components = entry.conjunction_of; + if ischar(components), components = {components}; end + A_all = []; + b_all = []; + for i = 1:numel(components) + comp = components{i}; + A_all = [A_all; pred.safety_limits.(comp).A_poly]; %#ok + b_all = [b_all; pred.safety_limits.(comp).b_poly]; %#ok + end + pred.mode_invariants.(name).A_poly = A_all; + pred.mode_invariants.(name).b_poly = b_all; + pred.mode_invariants.(name).meaning = entry.meaning; + pred.mode_invariants.(name).components = components; + end +end + +function group_out = parse_group(group_in, T_c0, T_f0, T_cold0, T_standby) + names = fieldnames(group_in); + group_out = struct(); for k = 1:numel(names) name = names{k}; - entry = J.predicates.(name); + % MATLAB jsondecode renames "_comment" -> "x_comment" and similar + if startsWith(name, '_') || startsWith(name, 'x_'), continue, end + entry = group_in.(name); + if ~isstruct(entry) || ~isfield(entry, 'halfspaces'), continue, end hs_list = entry.halfspaces; - - % jsondecode returns a struct array if entries are uniform, else cell. if iscell(hs_list) n_hs = numel(hs_list); get_hs = @(i) hs_list{i}; @@ -57,17 +89,19 @@ function pred = load_predicates(plant) n_hs = numel(hs_list); get_hs = @(i) hs_list(i); end - A_poly = zeros(n_hs, 10); b_poly = zeros(n_hs, 1); for i = 1:n_hs hs = get_hs(i); A_poly(i, hs.state_index) = hs.coeff; - b_poly(i) = eval(hs.rhs_expr); %#ok + b_poly(i) = evalin_context(hs.rhs_expr, T_c0, T_f0, T_cold0, T_standby); end - - pred.(name).A_poly = A_poly; - pred.(name).b_poly = b_poly; - pred.(name).meaning = entry.meaning; + group_out.(name).A_poly = A_poly; + group_out.(name).b_poly = b_poly; + group_out.(name).meaning = entry.meaning; end end + +function val = evalin_context(expr, T_c0, T_f0, T_cold0, T_standby) %#ok + val = eval(expr); +end diff --git a/reachability/predicates.json b/reachability/predicates.json index 8c5f0d9..30e46b8 100644 --- a/reachability/predicates.json +++ b/reachability/predicates.json @@ -1,24 +1,26 @@ { "_comment": [ - "Concretization of the FRET-spec predicates as numerical halfspaces.", - "This file is the single source of truth — all reach-analysis code loads", - "from here, and any future predicate changes happen here first.", + "Concretization of the FRET-spec predicates AND the hard safety limits.", + "Two categories kept distinct:", + " - operational_deadbands: soft bands around setpoint used by the DRC for", + " mode transitions (t_avg_in_range etc.). Violating these does not", + " cause damage, it just triggers a mode change or operator action.", + " - safety_limits: hard one-sided halfspaces corresponding to physical", + " damage mechanisms or reactor-trip setpoints. Barrier certificates", + " and reach-set safety checks should target THESE, not the deadbands.", "", - "Each predicate in fret-pipeline/pwr_hybrid_3.json (t_avg_above_min,", - "t_avg_in_range, p_above_crit, inv1_holds, inv2_holds) lives here as a", - "polytope over the 10-state vector x = [n, C1..C6, T_f, T_c, T_cold].", - "A polytope is {x : A_poly * x <= b_poly}, with units noted per entry." + "The FRET-spec invariants inv1_holds, inv2_holds are concretized as the", + "conjunction of relevant safety_limits for each mode." ], "_units": { "temperatures": "degrees Celsius (SI, internal model units)", - "n": "normalized power, 1.0 = full power", - "_display_note": "Figures and prints show Fahrenheit; predicates stored in C." + "n": "normalized power, 1.0 = full power" }, "references": { - "T_c0": "308.35 °C — full-power avg coolant (from pke_params.m)", - "T_f0": "328.35 °C — full-power fuel", - "T_cold0": "290.0 °C — full-power cold leg", - "T_standby": "275.02 °C — hot standby T_avg, defined as T_c0 - 33.33 C (= T_c0 - 60 F)" + "T_c0": "308.35 C — full-power avg coolant (from pke_params.m)", + "T_f0": "328.35 C — full-power fuel", + "T_cold0": "290.0 C — full-power cold leg", + "T_standby": "275.02 C — hot standby T_avg, defined as T_c0 - 33.33 C (= T_c0 - 60 F)" }, "derived": { "T_standby_offset_F": -60.0, @@ -27,63 +29,116 @@ "t_avg_in_range_halfwidth_C": 2.777777, "t_avg_above_min_margin_F": 10.0, "t_avg_above_min_margin_C": 5.555555, - "p_above_crit_threshold_n": 1.0e-4 + "p_above_crit_threshold_n": 1.0e-4, + "T_fuel_limit_C": 1200.0, + "T_c_high_trip_C": 320.0, + "n_high_trip": 1.15, + "T_cold_subcooling_margin_C": 15.0 }, - "predicates": { + + "operational_deadbands": { + "_comment": "Soft bands. Used by the DRC for mode switching, not for safety proofs.", "t_avg_above_min": { - "meaning": "Coolant has been warmed above a cold-start threshold; shutdown may transition to heatup.", - "concretization": "T_c >= T_standby + 5.556 C (= hot-standby + 10 F margin)", + "meaning": "Coolant warmed above cold-start threshold — shutdown may transition to heatup.", + "concretization": "T_c >= T_standby + 5.556 C (hot-standby + 10 F buffer)", "halfspaces": [ - { "state_index": 9, "coeff": -1.0, "rhs_expr": "-(T_standby + 5.556)", - "comment": "-T_c <= -(T_standby+5.556) i.e. T_c >= T_standby+5.556" } + { "state_index": 9, "coeff": -1.0, "rhs_expr": "-(T_standby + 5.556)" } ] }, "t_avg_in_range": { - "meaning": "Average coolant temperature is inside the operating band.", - "concretization": "|T_c - T_c0| <= 2.778 C (= 5 F tech-spec-like deadband)", + "meaning": "Average coolant in tight operating band — used for heatup->operation transition.", + "concretization": "|T_c - T_c0| <= 2.778 C (~5 F deadband)", "halfspaces": [ { "state_index": 9, "coeff": 1.0, "rhs_expr": "T_c0 + 2.778" }, { "state_index": 9, "coeff": -1.0, "rhs_expr": "-(T_c0 - 2.778)" } ] }, "p_above_crit": { - "meaning": "Reactor power is in the 'power range' instrumentation regime, above source/intermediate.", - "concretization": "n >= 1e-4 (0.01% of rated)", + "meaning": "Reactor power in the 'power range' instrumentation regime.", + "concretization": "n >= 1e-4", "halfspaces": [ - { "state_index": 1, "coeff": -1.0, "rhs_expr": "-1.0e-4", - "comment": "n >= 1e-4" } + { "state_index": 1, "coeff": -1.0, "rhs_expr": "-1.0e-4" } ] - }, - "inv1_holds": { - "meaning": "Heatup safety invariant — fuel not overheated, coolant subcooled, ramp rate nominal.", - "concretization": "(T_f <= 1200 C) AND (T_cold >= T_standby - 10 C) AND (dT_avg/dt <= 28 C/hr placeholder)", - "halfspaces": [ - { "state_index": 8, "coeff": 1.0, "rhs_expr": "1200.0", - "comment": "T_f <= 1200 C (fuel design limit; placeholder)" }, - { "state_index": 10, "coeff": -1.0, "rhs_expr": "-(T_standby - 10.0)", - "comment": "T_cold >= T_standby - 10 (minimal subcooling proxy)" } - ], - "_status": "PLACEHOLDER — ramp-rate constraint not expressible as a state halfspace without augmentation; DNBR not modeled." - }, - "inv2_holds": { - "meaning": "Operation safety invariant — power range, T_avg band, coolant subcooling.", - "concretization": "(n in [0.2, 1.1]) AND (T_c in [T_c0 - 8.33, T_c0 + 8.33] C) AND (T_cold in [T_cold0 - 20, T_cold0 + 10] C)", - "halfspaces": [ - { "state_index": 1, "coeff": 1.0, "rhs_expr": "1.1", "comment": "n <= 1.1" }, - { "state_index": 1, "coeff": -1.0, "rhs_expr": "-0.2", "comment": "n >= 0.2" }, - { "state_index": 9, "coeff": 1.0, "rhs_expr": "T_c0 + 8.33", "comment": "T_c <= T_c0 + 15 F" }, - { "state_index": 9, "coeff": -1.0, "rhs_expr": "-(T_c0 - 8.33)", "comment": "T_c >= T_c0 - 15 F" }, - { "state_index": 10, "coeff": 1.0, "rhs_expr": "T_cold0 + 10", "comment": "T_cold <= T_cold0 + 10" }, - { "state_index": 10, "coeff": -1.0, "rhs_expr": "-(T_cold0 - 20)", "comment": "T_cold >= T_cold0 - 20" } - ], - "_status": "PLACEHOLDER — DNBR not modeled; bands chosen to match typical PWR tech-spec deadbands but not calibrated against a specific plant." } }, + + "safety_limits": { + "_comment": [ + "Hard one-sided halfspaces. Exceeding any of these is damage or trip.", + "All are asymmetric — the plant is not equally vulnerable on both sides", + "of the setpoint. Values are representative of a 2-loop Westinghouse-", + "class PWR; calibrate to specific plant tech specs before defense." + ], + "fuel_centerline": { + "meaning": "Fuel centerline temperature below design limit to prevent UO2 melt.", + "concretization": "T_f <= 1200 C (conservative; actual melt ~2800 C)", + "halfspaces": [ + { "state_index": 8, "coeff": 1.0, "rhs_expr": "1200.0" } + ] + }, + "t_avg_high_trip": { + "meaning": "High-T_avg reactor trip. Typical PWR: ~612-616 F = 322-324 C.", + "concretization": "T_c <= 320 C (conservative)", + "halfspaces": [ + { "state_index": 9, "coeff": 1.0, "rhs_expr": "320.0" } + ] + }, + "t_avg_low_trip": { + "meaning": "Low-T_avg reactor trip. Typical PWR: ~540 F = 282 C.", + "concretization": "T_c >= 280 C", + "halfspaces": [ + { "state_index": 9, "coeff": -1.0, "rhs_expr": "-280.0" } + ] + }, + "n_high_trip": { + "meaning": "High-flux reactor trip. Typical PWR: 118%% of rated.", + "concretization": "n <= 1.15", + "halfspaces": [ + { "state_index": 1, "coeff": 1.0, "rhs_expr": "1.15" } + ] + }, + "n_low_operation": { + "meaning": "Operation mode is only valid at power (avoids intermediate-range instrumentation).", + "concretization": "n >= 0.15 (15%% of rated)", + "halfspaces": [ + { "state_index": 1, "coeff": -1.0, "rhs_expr": "-0.15" } + ] + }, + "cold_leg_subcooled": { + "meaning": "Cold leg stays subcooled with margin against loss-of-pressure events.", + "concretization": "T_cold <= T_cold0 + 15 (roughly saturation margin at operating pressure)", + "halfspaces": [ + { "state_index": 10, "coeff": 1.0, "rhs_expr": "T_cold0 + 15.0" } + ] + } + }, + + "mode_invariants": { + "_comment": [ + "Per-DRC-mode invariants: conjunctions of relevant safety_limits.", + "This is the target of per-mode reach and barrier analysis." + ], + "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." + }, + "inv2_holds": { + "meaning": "Operation mode safety envelope.", + "conjunction_of": [ + "fuel_centerline", + "t_avg_high_trip", + "t_avg_low_trip", + "n_high_trip", + "n_low_operation", + "cold_leg_subcooled" + ] + } + }, + "_placeholder_warning": [ - "The halfspace numbers for inv1_holds and inv2_holds are engineering", - "placeholders, not derived from a specific plant's tech specs or DNBR", - "correlation. They are the authors' best guess at the shape such", - "invariants should take and should be revised before the thesis defense.", - "The t_avg_* and p_above_crit predicates are defensible for the demo." + "Numerical values in safety_limits are representative (2-loop Westinghouse-", + "class PWR tech-spec ranges) but NOT calibrated to a specific plant.", + "Calibrate against a real plant's tech specs before defense." ] } diff --git a/reachability/reach_operation.m b/reachability/reach_operation.m index 39d1f33..0dd854f 100644 --- a/reachability/reach_operation.m +++ b/reachability/reach_operation.m @@ -81,9 +81,13 @@ Q_max = 1.00 * plant.P0; dQ_lo = Q_min - Q_nom; % -0.15 * P0 dQ_hi = Q_max - Q_nom; % 0 -% X_safe: pulled directly from the t_avg_in_range predicate (the -% reachability/predicates.json concretization). Row 1 of its A_poly is -% [0,...,0,+1,0] with b = T_c0 + halfwidth, row 2 is the lower bound. +% X_safe is inv2_holds — the operation-mode safety envelope. Each row of +% inv2.A_poly is a hard safety limit (fuel centerline, T_c high/low trip, +% n high/low, cold-leg subcooling). We check reach-tube containment +% halfspace-by-halfspace so failure modes are attributable. +inv2 = pred.mode_invariants.inv2_holds; +% Keep the old +-halfwidth report too, for continuity with the +% operational-deadband framing. delta_safe_Tc = pred.constants.t_avg_in_range_halfwidth_C; % [C] %% ===== Reach set ===== @@ -119,6 +123,26 @@ else fprintf(' OK: reach set stays inside the safe band.\n'); end +%% Hard safety-limit check (inv2_holds halfspace-by-halfspace) +% For each row a_k of inv2.A_poly with threshold b_k, check whether +% max over reach tube of a_k * x stays <= b_k. The reach tube upper +% envelope is Xabs_hi; lower envelope is Xabs_lo. We evaluate +% max(a_k * x) using Xabs_hi where a_k > 0, Xabs_lo where a_k < 0. +fprintf('\n=== Operation-mode reach vs inv2_holds safety limits ===\n'); +A_inv = inv2.A_poly; b_inv = inv2.b_poly; +comps = inv2.components; +for k = 1:size(A_inv, 1) + a = A_inv(k, :).'; + % envelope maximum of a' * x across the reach tube + x_envelope = Xabs_hi .* (a > 0) + Xabs_lo .* (a < 0); % 10 x M + max_ax = max(a.' * x_envelope); + margin = b_inv(k) - max_ax; + status = 'OK'; + if margin < 0, status = '*** VIOLATED ***'; end + fprintf(' [%s] a''x <= %.3f | max a''x = %.3f | margin = %+.3f %s\n', ... + comps{k}, b_inv(k), max_ax, margin, status); +end + %% Per-state reach-set growth diagnostic (final time vs initial) state_names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'}; fprintf('\n=== Reach-set width at t=0 vs t=T_final ===\n');