danesabo/PWR-HYBRID-3
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

predicates: split operational deadbands from hard safety limits

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
Dane Sabo 2026-04-20 16:04:40 -04:00
parent e69fd0a6f4
commit a20d2a05e9
6 changed files with 233 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
docs/figures/reach_operation_Tc.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 62 KiB

After

Width:  |  Height:  |  Size: 62 KiB

BIN
docs/figures/reach_operation_n.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 34 KiB

After

Width:  |  Height:  |  Size: 34 KiB

34
reachability/barrier_lyapunov.m
View File

@ -52,6 +52,22 @@ P = lyap(A_cl.', Qbar);
assert(all(eig(P) > 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');

86
reachability/load_predicates.m
View File

@ -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'));
@ -40,16 +36,52 @@ function pred = load_predicates(plant)
'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<AGROW>
b_all = [b_all; pred.safety_limits.(comp).b_poly]; %#ok<AGROW>
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<EVLDF>
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<INUSD>
val = eval(expr);
end

163
reachability/predicates.json
View File

@ -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."
]
}

30
reachability/reach_operation.m
View File

@ -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');