diff --git a/docs/figures/linearize_sanity.png b/docs/figures/linearize_sanity.png index d5528fe..a4b3c93 100644 Binary files a/docs/figures/linearize_sanity.png and b/docs/figures/linearize_sanity.png differ diff --git a/docs/figures/mode_sweep_1_shutdown.png b/docs/figures/mode_sweep_1_shutdown.png index 1068d7f..bc75195 100644 Binary files a/docs/figures/mode_sweep_1_shutdown.png and b/docs/figures/mode_sweep_1_shutdown.png differ diff --git a/docs/figures/mode_sweep_2_heatup.png b/docs/figures/mode_sweep_2_heatup.png index 7de1db4..9697c19 100644 Binary files a/docs/figures/mode_sweep_2_heatup.png and b/docs/figures/mode_sweep_2_heatup.png differ diff --git a/docs/figures/mode_sweep_3_operation.png b/docs/figures/mode_sweep_3_operation.png index aef2535..7f401ab 100644 Binary files a/docs/figures/mode_sweep_3_operation.png and b/docs/figures/mode_sweep_3_operation.png differ diff --git a/docs/figures/mode_sweep_3b_operation_lqr.png b/docs/figures/mode_sweep_3b_operation_lqr.png index 7907679..aed81a6 100644 Binary files a/docs/figures/mode_sweep_3b_operation_lqr.png and b/docs/figures/mode_sweep_3b_operation_lqr.png differ diff --git a/docs/figures/mode_sweep_4_scram.png b/docs/figures/mode_sweep_4_scram.png index a7c8bea..219b952 100644 Binary files a/docs/figures/mode_sweep_4_scram.png and b/docs/figures/mode_sweep_4_scram.png differ diff --git a/docs/figures/mode_sweep_heatup_tracking.png b/docs/figures/mode_sweep_heatup_tracking.png index 0e9530e..1323c9b 100644 Binary files a/docs/figures/mode_sweep_heatup_tracking.png and b/docs/figures/mode_sweep_heatup_tracking.png differ diff --git a/docs/figures/mode_sweep_op_P_vs_LQR.png b/docs/figures/mode_sweep_op_P_vs_LQR.png index 046eb86..79d1c69 100644 Binary files a/docs/figures/mode_sweep_op_P_vs_LQR.png and b/docs/figures/mode_sweep_op_P_vs_LQR.png differ diff --git a/docs/figures/mode_sweep_power_overview.png b/docs/figures/mode_sweep_power_overview.png index 479e6fa..89f6734 100644 Binary files a/docs/figures/mode_sweep_power_overview.png and b/docs/figures/mode_sweep_power_overview.png differ diff --git a/docs/figures/reach_operation_Tc.png b/docs/figures/reach_operation_Tc.png index d918ba3..ee90c36 100644 Binary files a/docs/figures/reach_operation_Tc.png and b/docs/figures/reach_operation_Tc.png differ diff --git a/docs/figures/reach_operation_n.png b/docs/figures/reach_operation_n.png index 44d1e11..07d3981 100644 Binary files a/docs/figures/reach_operation_n.png and b/docs/figures/reach_operation_n.png differ diff --git a/plant-model/main_mode_sweep.m b/plant-model/main_mode_sweep.m index 95deef8..bbe33c0 100644 --- a/plant-model/main_mode_sweep.m +++ b/plant-model/main_mode_sweep.m @@ -15,8 +15,10 @@ clear; clc; close all; addpath('controllers'); +addpath(fullfile('..', 'reachability')); % for load_predicates plant = pke_params(); +pred = load_predicates(plant); % T_standby and FRET-predicate concretizations figdir = fullfile('..', 'docs', 'figures'); if ~exist(figdir, 'dir'), mkdir(figdir); end @@ -26,24 +28,23 @@ CtoF = @(T) T*9/5 + 32; % Operating IC: full-power steady state. x0_op = pke_initial_conditions(plant); -% Hot-standby shutdown IC: coolant everywhere at T_cold0 (290 C), reactor -% deep subcritical with only a trace neutron population. Note: our -% feedback model is referenced to the hot full-power state, so going -% below ~290 C violates the linear-coefficient assumption. This IC is -% the lowest T we trust the model at. +% Hot-standby shutdown IC: coolant flat at T_standby = T_c0 - 60 F ~ 275 C, +% reactor deep subcritical with only a trace neutron population. T_standby +% comes from reachability/predicates.json (single source of truth). Well +% inside the model's +/-50 C trust region around the operating point, and +% above coolant saturation at reduced operating pressure. +T_standby = pred.constants.T_standby; n_shut = 1e-6; C_shut = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_shut; -x0_shut = [n_shut; C_shut; plant.T_cold0; plant.T_cold0; plant.T_cold0]; +x0_shut = [n_shut; C_shut; T_standby; T_standby; T_standby]; -% Heatup IC: reactor already taken critical at 0.1% power, low-power -% criticality achieved in shutdown mode before the DRC transitions. -% This avoids the unphysical "ramp from n=0" scenario where the P -% controller has to build power from decay heat before temperature can -% move at all. Mirrors real plant practice: achieve criticality, then -% heat up. +% Heatup IC: reactor already taken critical at 0.1% power at hot-standby +% temperature. Mirrors real plant practice: achieve criticality, then +% heat up. Same T_standby as the shutdown IC — heatup begins from where +% shutdown left off. n_heat = 1e-3; C_heat = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_heat; -x0_heat = [n_heat; C_heat; plant.T_cold0; plant.T_cold0; plant.T_cold0]; +x0_heat = [n_heat; C_heat; T_standby; T_standby; T_standby]; %% ===== Mode 1: SHUTDOWN ===== fprintf('\n===== Mode 1: ctrl_shutdown =====\n'); @@ -54,11 +55,11 @@ tspan_shut = [0, 600]; %% ===== Mode 2: HEATUP ===== fprintf('\n===== Mode 2: ctrl_heatup =====\n'); ref_heatup = struct(); -ref_heatup.T_start = plant.T_cold0; % 290 C -ref_heatup.T_target = plant.T_c0; % 308.35 C +ref_heatup.T_start = T_standby; % ~275 C (hot standby, = T_c0 - 60 F) +ref_heatup.T_target = plant.T_c0; % 308.35 C (operating setpoint) ref_heatup.ramp_rate = 28 / 3600; % 28 C/hr, tech-spec limit Q_sg_heat = @(t) 0; % no SG demand during heatup -tspan_heat = [0, 3000]; % ~50 min +tspan_heat = [0, 5400]; % ~90 min (33 C span at 28 C/hr = 71 min + settling) [t2, X2, U2] = pke_solver(plant, Q_sg_heat, @ctrl_heatup, ref_heatup, tspan_heat, x0_heat); %% ===== Mode 3a: OPERATION (plain P) ===== diff --git a/reachability/load_predicates.m b/reachability/load_predicates.m new file mode 100644 index 0000000..79bc6b0 --- /dev/null +++ b/reachability/load_predicates.m @@ -0,0 +1,73 @@ +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 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. +% +% 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 + + here = fileparts(mfilename('fullpath')); + raw = fileread(fullfile(here, 'predicates.json')); + J = jsondecode(raw); + + % --- Constants used in rhs_expr evaluations --- + T_c0 = plant.T_c0; %#ok + T_f0 = plant.T_f0; %#ok + T_cold0 = plant.T_cold0; %#ok + T_standby_offset_C = J.derived.T_standby_offset_C; + 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_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); + + % --- Loop over predicates, build A/b matrices --- + names = fieldnames(J.predicates); + for k = 1:numel(names) + name = names{k}; + entry = J.predicates.(name); + 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}; + else + 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 + end + + pred.(name).A_poly = A_poly; + pred.(name).b_poly = b_poly; + pred.(name).meaning = entry.meaning; + end +end diff --git a/reachability/predicates.json b/reachability/predicates.json new file mode 100644 index 0000000..8c5f0d9 --- /dev/null +++ b/reachability/predicates.json @@ -0,0 +1,89 @@ +{ + "_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.", + "", + "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." + ], + "_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." + }, + "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)" + }, + "derived": { + "T_standby_offset_F": -60.0, + "T_standby_offset_C": -33.333333333, + "t_avg_in_range_halfwidth_F": 5.0, + "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 + }, + "predicates": { + "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)", + "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" } + ] + }, + "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)", + "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)", + "halfspaces": [ + { "state_index": 1, "coeff": -1.0, "rhs_expr": "-1.0e-4", + "comment": "n >= 1e-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." + } + }, + "_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." + ] +} diff --git a/reachability/reach_operation.m b/reachability/reach_operation.m index 183d046..39d1f33 100644 --- a/reachability/reach_operation.m +++ b/reachability/reach_operation.m @@ -40,6 +40,7 @@ addpath('../plant-model', '../plant-model/controllers'); plant = pke_params(); x_op = pke_initial_conditions(plant); +pred = load_predicates(plant); % single source of truth for predicate bands %% ===== Closed-loop linearization ===== [A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); @@ -80,9 +81,10 @@ Q_max = 1.00 * plant.P0; dQ_lo = Q_min - Q_nom; % -0.15 * P0 dQ_hi = Q_max - Q_nom; % 0 -% X_safe: |T_c - T_c0| <= 5 K (roughly matching the "t_avg_in_range" -% predicate window — adjust once the FRET predicate threshold is locked). -delta_safe_Tc = 5.0; % [C] +% 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. +delta_safe_Tc = pred.constants.t_avg_in_range_halfwidth_C; % [C] %% ===== Reach set ===== % Propagate in deviation coordinates: dx = x - x_op.