From d2997c2861684133e74a3f9813794e92547fd052 Mon Sep 17 00:00:00 2001 From: Dane Sabo Date: Fri, 17 Apr 2026 12:52:03 -0400 Subject: [PATCH] plant-model: add shutdown/heatup/scram controllers and LQR, linearize MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fill out the DRC mode set with ctrl_shutdown (u = -5*beta), ctrl_scram (u = -8*beta), and ctrl_heatup (feedback-linearizing P on ramped T_avg reference, saturated u, no integrator). Add ctrl_operation_lqr as a full-state-feedback counterpart to ctrl_operation — K cached, closed-loop essentially perfect under the 100%->80% Q_sg step where plain P has ~5F overshoot. Add pke_linearize for numerical (A, B, B_w) Jacobians at any operating point; test_linearize confirms ~4e-4 rel err vs nonlinear sim for a 5% Q_sg step. Extend pke_solver with an optional x0 argument so each mode can start from a plausible IC. main_mode_sweep.m exercises all five modes back-to-back and saves the 4-panel plots. CLAUDE.md updated with model-validity-range note (trust region is ~+/-50C around operating point; true cold shutdown is out of scope for the linear feedback coefficients). Hacker-Split: build out control layer end-to-end for reachability. Co-Authored-By: Claude Opus 4.7 (1M context) --- plant-model/CLAUDE.md | 56 +++++++- plant-model/controllers/ctrl_heatup.m | 57 ++++++++ plant-model/controllers/ctrl_operation_lqr.m | 53 +++++++ plant-model/controllers/ctrl_scram.m | 20 +++ plant-model/controllers/ctrl_shutdown.m | 21 +++ plant-model/main_mode_sweep.m | 137 +++++++++++++++++++ plant-model/pke_linearize.m | 54 ++++++++ plant-model/pke_solver.m | 8 +- plant-model/test_linearize.m | 106 ++++++++++++++ 9 files changed, 503 insertions(+), 9 deletions(-) create mode 100644 plant-model/controllers/ctrl_heatup.m create mode 100644 plant-model/controllers/ctrl_operation_lqr.m create mode 100644 plant-model/controllers/ctrl_scram.m create mode 100644 plant-model/controllers/ctrl_shutdown.m create mode 100644 plant-model/main_mode_sweep.m create mode 100644 plant-model/pke_linearize.m create mode 100644 plant-model/test_linearize.m diff --git a/plant-model/CLAUDE.md b/plant-model/CLAUDE.md index f8665b4..b7fdf30 100644 --- a/plant-model/CLAUDE.md +++ b/plant-model/CLAUDE.md @@ -62,14 +62,20 @@ plant-model/ pke_params.m plant parameters and derived steady state pke_initial_conditions.m analytic equilibrium x0 from plant pke_th_rhs.m dynamics: dx/dt = f(t, x, plant, Q_sg, u) - pke_solver.m closed-loop driver; takes a controller fn handle + pke_linearize.m numerical (A, B, B_w) at any x_star + pke_solver.m closed-loop driver (optional x0 arg) plot_pke_results.m 4-panel results plot load_profile.m canned Q_sg demand shapes - main.m entry point — scenario config + run + main.m entry point — single-mode scenarios + main_mode_sweep.m runs all DRC modes back-to-back + test_linearize.m sanity-check Jacobians vs nonlinear sim controllers/ ctrl_null.m u = 0 (baseline, no feedback) + ctrl_shutdown.m hot standby: u = -5*beta (constant) + ctrl_heatup.m ramp T_avg: feedback lin. + P on ramp, saturated u ctrl_operation.m stabilizing: P on T_avg - (future: ctrl_heatup, ctrl_scram, ctrl_shutdown) + ctrl_operation_lqr.m full-state LQR around x_op (persistent-cached K) + ctrl_scram.m emergency: u = -8*beta (constant) ``` **Controller contract:** @@ -117,16 +123,52 @@ interior solver evaluations — revisit if / when that matters. - Controllers are pure functions — no persistent state. If a mode ever needs integral action, augment `x` rather than hiding state in globals. -## For reachability work (next chapter) +## Model validity range + +The feedback coefficients `alpha_f`, `alpha_c` are *linear* slopes taken +about the full-power operating point `(T_f0, T_c0)`. They become +unphysical when extrapolated far from that reference — e.g. at true +cold shutdown (~50 C everywhere), the model predicts ~+5 $ intrinsic +reactivity that real reactors do not exhibit because temperature +coefficients saturate and flip sign in the cold-unborated state. + +**Trust range:** roughly ±50 C around the operating point. "Shutdown" +in this demo means hot standby (~290 C, n≈0), not cold shutdown. +Extending to cold-shutdown semantics requires piecewise-linear or +table-lookup temperature coefficients and is out of scope for the +running example. + +## For reachability work - `pke_th_rhs(t, x, plant, Q_sg, u)` is the `ẋ = f(x, u)` with `Q_sg` as the disturbance `w`. For a given mode, substitute `u = ctrl_(t, x, plant, ref)` to get the closed-loop `ẋ = f_cl(x, w)`. -- Linearize around `x0` from `pke_initial_conditions(plant)` for `A, B, B_w` - matrices when a linear reach set is enough. +- Use `pke_linearize(plant, x_star, u_star, Q_star)` for `(A, B, B_w)` at + any operating point; returns central-finite-difference Jacobians. The + result at the full-power point is saved to + `../reachability/linearization_at_op.mat` by `test_linearize.m`. - Safe regions come from the FRET spec's guards (see `../fret-pipeline/specs/synthesis_config_v3.json`) — these define - `X_entry`, `X_exit`, `X_safe` per the thesis approach. + `X_entry`, `X_exit`, `X_safe` per the thesis approach. Numerical + predicate thresholds still need to be pinned (see + `../claude_memory/` notes for the discussion). - Do NOT try to verify the whole hybrid system at once. Pick one mode, compute its reach set, discharge the per-mode obligation. The compositionality argument in the thesis is what makes this tractable. +- First reach set is in `../reachability/reach_operation.m` (operation + mode under LQR); a Lyapunov-ellipsoid barrier-cert attempt is in + `../reachability/barrier_lyapunov.m`. + +## For control design + +- The LQR gain in `ctrl_operation_lqr.m` is cached in a persistent + variable. If you change `Q_lqr` or `R_lqr` or mutate `plant`, restart + MATLAB to clear the cache. +- `ctrl_heatup.m` uses feedback linearization + saturated P. No + integrator on purpose — adding PI would double the plant's intrinsic + thermal-mass integration and make anti-windup a hybrid mode, both + bad for reach. +- Starting heatup from a truly zero-power IC produces a large lag + + power-spike overshoot in the simulation. Physical resolution: the + reactor should be taken critical in shutdown mode (at ~0.1% power) + before DRC transitions to heatup. `main_mode_sweep.m` uses this IC. diff --git a/plant-model/controllers/ctrl_heatup.m b/plant-model/controllers/ctrl_heatup.m new file mode 100644 index 0000000..eca79a8 --- /dev/null +++ b/plant-model/controllers/ctrl_heatup.m @@ -0,0 +1,57 @@ +function u = ctrl_heatup(t, x, plant, ref) +% CTRL_HEATUP Ramp T_avg toward a target at a bounded rate. +% +% Structure: +% u_ff = -alpha_f*(T_f - T_f0) - alpha_c*(T_c - T_c0) % cancel feedback +% T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target) % ramped reference +% u_unsat = u_ff + Kp * (T_ref - T_avg) % P on error +% u = sat(u_unsat, ref.u_min, ref.u_max) % bounded rod worth +% +% Why saturation: +% Without it the P gain can push u toward prompt-supercritical as the +% cold-hot feedback bias unwinds late in the ramp. Capping u at +% +0.5*beta guarantees rho_total < beta (below prompt), which in +% turn bounds the neutron-kinetics excursion rate for reachability. +% +% Why no integrator: +% Ramp tracking has a structural lag proportional to ramp_rate / Kp_eff. +% Acceptable because the DRC exits heatup on a predicate window +% (t_avg_in_range & p_above_crit), not on zero steady-state error. +% Adding PI would double-count the intrinsic plant integrator +% (thermal mass) and make anti-windup a hybrid transition. +% +% Inputs: +% t - time [s] +% x - state vector (10 x 1) +% plant - parameter struct (alpha_f, alpha_c, T_f0, T_c0 used) +% ref - struct with fields: +% .T_start starting T_avg [C] +% .T_target final T_avg [C] +% .ramp_rate desired dT_avg/dt [C/s] +% .u_min (optional) lower saturation [dk/k]; default -5*beta +% .u_max (optional) upper saturation [dk/k]; default +0.5*beta + + Kp = 1e-4; % [dk/k per K] + + T_f = x(8); + T_c = x(9); + T_avg = T_c; + + % Feedforward: cancel intrinsic temperature feedback so rho_total = Kp*e + % (before saturation). + u_ff = -plant.alpha_f * (T_f - plant.T_f0) ... + -plant.alpha_c * (T_avg - plant.T_c0); + + % Ramped reference, clamped at target. + T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target); + + e = T_ref - T_avg; + u_unsat = u_ff + Kp * e; + + % Saturation bounds (defaults keep rod worth subcritical-prompt). + if isfield(ref, 'u_min'), u_min = ref.u_min; else, u_min = -5 * plant.beta; end + if isfield(ref, 'u_max'), u_max = ref.u_max; else, u_max = 0.5 * plant.beta; end + + u = min(max(u_unsat, u_min), u_max); + +end diff --git a/plant-model/controllers/ctrl_operation_lqr.m b/plant-model/controllers/ctrl_operation_lqr.m new file mode 100644 index 0000000..cd90691 --- /dev/null +++ b/plant-model/controllers/ctrl_operation_lqr.m @@ -0,0 +1,53 @@ +function u = ctrl_operation_lqr(t, x, plant, ref) +% CTRL_OPERATION_LQR Full-state feedback for operation mode. +% +% u = -K * (x - x_op), K from LQR on (A, B) at x_op. +% +% Advantages over ctrl_operation (plain P on T_avg): +% - Damps the fast neutronics and the slow thermal response jointly. +% - Adds NO state (K is 1x10, feedback is memoryless). +% - Closed-loop A_cl = A - B*K is trivially propagable for reachability. +% +% Q/R are persistently cached on the first call. If the caller mutates +% plant after the first invocation, gains become stale — clear `persistent` +% by restarting MATLAB. That's a reach-analysis bookkeeping issue, not a +% runtime one. +% +% Weighting rationale: +% Q diagonal: +% n -> 1 (power tracking matters) +% C1..C6 -> 1e-3 (precursor deviations are informational) +% T_f -> 1e-2 (care about it but not as much as coolant) +% T_c -> 1e2 (primary objective) +% T_cold -> 1 (secondary but consequential) +% R = 1e6 (discourage large rod movements — one Kp of |u| ~ 1e-3) +% +% Inputs: +% t, x, plant, ref - standard signature. ref is unused; the setpoint +% is baked in as x_op. If you want to track a +% different operating point, regenerate K there. + + persistent K x_op_cached + + if isempty(K) + x_op_cached = pke_initial_conditions(plant); + [A, B, ~] = pke_linearize(plant, x_op_cached, 0, plant.P0); + + Q = diag([1, ... + 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, ... + 1e-2, 1e2, 1]); + R = 1e6; + + % MATLAB's lqr requires Control System Toolbox. Fall back to + % icare (in base MATLAB from R2019a) if lqr is unavailable. + try + K = lqr(A, B, Q, R); + catch + [~, ~, K_ric] = icare(A, B, Q, R); + K = K_ric; + end + end + + u = -K * (x - x_op_cached); + +end diff --git a/plant-model/controllers/ctrl_scram.m b/plant-model/controllers/ctrl_scram.m new file mode 100644 index 0000000..b07b72a --- /dev/null +++ b/plant-model/controllers/ctrl_scram.m @@ -0,0 +1,20 @@ +function u = ctrl_scram(t, x, plant, ref) +% CTRL_SCRAM Emergency shutdown: max negative rod worth inserted. +% +% Applies the full scram rod worth instantaneously. In a real plant the +% rods free-fall in ~2-3 seconds; this idealized version steps to the +% final worth at t = t_scram. Good enough for first-pass reachability, +% conservative for safety arguments (faster insertion => less excursion). +% +% u = -8 * beta (~8 $, in the typical 8-10 $ regulatory band) +% +% Same feedback-linearization comment as ctrl_shutdown: this dominates +% any plausible feedback the reactor can produce en route from any +% power/temperature state. +% +% Inputs: +% t, x, plant, ref - standard signature; all unused here. + + u = -8 * plant.beta; + +end diff --git a/plant-model/controllers/ctrl_shutdown.m b/plant-model/controllers/ctrl_shutdown.m new file mode 100644 index 0000000..a56014b --- /dev/null +++ b/plant-model/controllers/ctrl_shutdown.m @@ -0,0 +1,21 @@ +function u = ctrl_shutdown(t, x, plant, ref) +% CTRL_SHUTDOWN Hot-standby / cold-shutdown passive controller. +% +% Holds the reactor deeply subcritical with rods fully in. No feedback — +% just a constant negative external reactivity that dominates any +% temperature-feedback contribution the reactor could plausibly produce +% in its accessible state set. +% +% u = -5 * beta (~5 $ subcritical) +% +% Sign check: rho_total = u + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0). +% At cold IC (T < T0) the feedback contribution is positive (roughly +% +2.8 $ at T = T_cold0 everywhere), so -5 $ leaves the reactor at about +% -2.2 $ — comfortably subcritical. +% +% Inputs: +% t, x, plant, ref - standard signature; all unused here. + + u = -5 * plant.beta; + +end diff --git a/plant-model/main_mode_sweep.m b/plant-model/main_mode_sweep.m new file mode 100644 index 0000000..95deef8 --- /dev/null +++ b/plant-model/main_mode_sweep.m @@ -0,0 +1,137 @@ +%% main_mode_sweep.m — demo each DRC continuous-mode controller +% +% One scenario per mode, from an initial condition plausible for that +% mode. Figures saved to ../docs/figures/mode_sweep_*.png so the runs +% are reviewable without the MATLAB IDE. +% +% Modes, in the order DRC visits them starting from shutdown: +% 1. shutdown — deep subcritical, cold IC, Q_sg = 0 +% 2. heatup — same cold IC, ramped T_avg reference +% 3. operation — operating IC, 100% -> 80% Q_sg step +% 4. scram — operating IC, rods slammed in, Q_sg decays to 0 +% +% Each run produces plot_pke_results 4-panel plus saves a lightweight +% summary figure. + +clear; clc; close all; +addpath('controllers'); + +plant = pke_params(); +figdir = fullfile('..', 'docs', 'figures'); +if ~exist(figdir, 'dir'), mkdir(figdir); end + +CtoF = @(T) T*9/5 + 32; + +%% ===== IC helpers ===== +% 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. +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]; + +% 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. +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]; + +%% ===== Mode 1: SHUTDOWN ===== +fprintf('\n===== Mode 1: ctrl_shutdown =====\n'); +Q_sg_shut = @(t) 0; % no SG demand +tspan_shut = [0, 600]; +[t1, X1, U1] = pke_solver(plant, Q_sg_shut, @ctrl_shutdown, [], tspan_shut, x0_shut); + +%% ===== 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.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 +[t2, X2, U2] = pke_solver(plant, Q_sg_heat, @ctrl_heatup, ref_heatup, tspan_heat, x0_heat); + +%% ===== Mode 3a: OPERATION (plain P) ===== +fprintf('\n===== Mode 3a: ctrl_operation (P) =====\n'); +Q_sg_op = @(t) plant.P0 * (1.0 - 0.2 * (t >= 30)); % 100% -> 80% step +ref_op = struct('T_avg', plant.T_c0); +tspan_op = [0, 600]; +[t3, X3, U3] = pke_solver(plant, Q_sg_op, @ctrl_operation, ref_op, tspan_op, x0_op); + +%% ===== Mode 3b: OPERATION (LQR) ===== +fprintf('\n===== Mode 3b: ctrl_operation_lqr =====\n'); +[t3b, X3b, U3b] = pke_solver(plant, Q_sg_op, @ctrl_operation_lqr, [], tspan_op, x0_op); + +%% ===== Mode 4: SCRAM ===== +fprintf('\n===== Mode 4: ctrl_scram =====\n'); +% After a scram signal, the turbine trips and SG isolation occurs fast. +% Q_sg drops from P0 to a decay-heat-level sink (3% of P0) in ~10 s, +% then holds. Not a true decay-heat model — good enough to show the +% post-scram thermal trajectory without coolant going unphysically cold. +Q_sg_scr = @(t) plant.P0 * max(0.03, 1 - max(0, t - 10) / 10); +tspan_scr = [0, 600]; +[t4, X4, U4] = pke_solver(plant, Q_sg_scr, @ctrl_scram, [], tspan_scr, x0_op); + +%% ===== Per-mode 4-panel plots ===== +plot_pke_results(t1, X1, U1, plant, Q_sg_shut, 'ctrl\_shutdown (cold IC)'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_1_shutdown.png'), 'Resolution', 150); + +plot_pke_results(t2, X2, U2, plant, Q_sg_heat, 'ctrl\_heatup (ramp T\_avg)'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_2_heatup.png'), 'Resolution', 150); + +plot_pke_results(t3, X3, U3, plant, Q_sg_op, 'ctrl\_operation (P on T\_avg)'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3_operation.png'), 'Resolution', 150); + +plot_pke_results(t3b, X3b, U3b, plant, Q_sg_op, 'ctrl\_operation\_lqr'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3b_operation_lqr.png'), 'Resolution', 150); + +plot_pke_results(t4, X4, U4, plant, Q_sg_scr, 'ctrl\_scram'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_4_scram.png'), 'Resolution', 150); + +%% ===== Heatup ramp-tracking figure ===== +% Overlay the T_ref signal on T_avg so the lag is visible. +T_ref_trace = min(ref_heatup.T_start + ref_heatup.ramp_rate .* t2, ref_heatup.T_target); +figure('Position', [100 80 900 400], 'Name', 'Heatup ramp tracking'); +plot(t2/60, CtoF(T_ref_trace), 'k--', 'LineWidth', 1.2); hold on; +plot(t2/60, CtoF(X2(:,9)), 'r-', 'LineWidth', 1.5); +xlabel('Time [min]'); ylabel('T_{avg} [F]'); +legend('T_{ref}(t) (28 C/hr ramp)', 'T_{avg} (plant)', 'Location', 'southeast'); +title('ctrl\_heatup: reference tracking'); +grid on; +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_heatup_tracking.png'), 'Resolution', 150); + +%% ===== P vs LQR head-to-head on T_avg ===== +figure('Position', [100 80 900 400], 'Name', 'Operation: P vs LQR'); +plot(t3, CtoF(X3(:,9)), 'b-', 'LineWidth', 1.5); hold on; +plot(t3b, CtoF(X3b(:,9)), 'r-', 'LineWidth', 1.5); +yline(CtoF(plant.T_c0), 'k--'); +xlabel('Time [s]'); ylabel('T_{avg} [F]'); +legend('ctrl\_operation (P)', 'ctrl\_operation\_lqr', 'setpoint', 'Location', 'best'); +title('Operation mode: P vs LQR under 100% \rightarrow 80% Q_{sg} step'); +grid on; +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_op_P_vs_LQR.png'), 'Resolution', 150); + +%% ===== Summary figure: normalized power across all modes ===== +figure('Position', [100 80 1000 450], 'Name', 'Normalized power, all modes'); +subplot(2,2,1); semilogy(t1, max(X1(:,1), 1e-20)); grid on; + title('Shutdown'); xlabel('t [s]'); ylabel('n (log)'); +subplot(2,2,2); plot(t2/60, X2(:,1)); grid on; + title('Heatup'); xlabel('t [min]'); ylabel('n'); +subplot(2,2,3); plot(t3, X3(:,1)); grid on; + title('Operation'); xlabel('t [s]'); ylabel('n'); +subplot(2,2,4); semilogy(t4, max(X4(:,1), 1e-20)); grid on; + title('Scram'); xlabel('t [s]'); ylabel('n (log)'); +sgtitle('Normalized reactor power n(t) across DRC modes'); +exportgraphics(gcf, fullfile(figdir, 'mode_sweep_power_overview.png'), 'Resolution', 150); + +fprintf('\n=== Figures written to %s ===\n', figdir); diff --git a/plant-model/pke_linearize.m b/plant-model/pke_linearize.m new file mode 100644 index 0000000..7eb60fc --- /dev/null +++ b/plant-model/pke_linearize.m @@ -0,0 +1,54 @@ +function [A, B, B_w, x_star, u_star, Q_star] = pke_linearize(plant, x_star, u_star, Q_star) +% PKE_LINEARIZE Numerical Jacobians of the PKE+T/H RHS at an operating point. +% +% Computes (A, B, B_w) such that, for small perturbations dx, du, dw +% around (x_star, u_star, Q_star): +% +% dx/dt ~ A*dx + B*du + B_w*dw, where w = Q_sg. +% +% Uses central finite differences. Step sizes are scaled to each +% variable's magnitude; fall back to an absolute epsilon near zero. +% +% Inputs: +% plant - parameter struct from pke_params() +% x_star - (optional) 10x1 linearization point; default = operating x0 +% u_star - (optional) scalar control; default = 0 +% Q_star - (optional) scalar SG heat removal [W]; default = plant.P0 +% +% Outputs: +% A - 10x10 state Jacobian df/dx +% B - 10x1 input Jacobian df/du +% B_w - 10x1 disturbance Jacobian df/d(Q_sg) +% x_star, u_star, Q_star - echoed back for the caller's records + + if nargin < 2 || isempty(x_star), x_star = pke_initial_conditions(plant); end + if nargin < 3 || isempty(u_star), u_star = 0; end + if nargin < 4 || isempty(Q_star), Q_star = plant.P0; end + + n = numel(x_star); + eps_rel = 1e-6; + eps_abs = 1e-8; + + % Wrap the RHS as a function of (x, u, w) with plant fixed. + f = @(x, u, w) pke_th_rhs(0, x, plant, @(t) w, u); + + f0 = f(x_star, u_star, Q_star); %#ok % sanity slot + + % --- A: df/dx via central differences, column by column --- + A = zeros(n, n); + for k = 1:n + h = max(eps_rel * abs(x_star(k)), eps_abs); + xp = x_star; xp(k) = xp(k) + h; + xm = x_star; xm(k) = xm(k) - h; + A(:, k) = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) / (2*h); + end + + % --- B: df/du --- + h = max(eps_rel * abs(u_star), eps_abs); + B = (f(x_star, u_star + h, Q_star) - f(x_star, u_star - h, Q_star)) / (2*h); + + % --- B_w: df/dQ_sg --- + h = max(eps_rel * abs(Q_star), 1.0); % Q_star ~ 1e9, abs floor at 1 W + B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) / (2*h); + +end diff --git a/plant-model/pke_solver.m b/plant-model/pke_solver.m index d499ad2..12291d7 100644 --- a/plant-model/pke_solver.m +++ b/plant-model/pke_solver.m @@ -1,4 +1,4 @@ -function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan) +function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan, x0) % PKE_SOLVER Solve the coupled PKE + T/H system in closed loop with a mode controller. % % Inputs: @@ -7,6 +7,8 @@ function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan) % ctrl_fn - function handle u = ctrl_fn(t, x, plant, ref) % ref - struct of per-mode setpoints (e.g. ref.T_avg); [] if unused % tspan - [t_start, t_end] in seconds +% x0 - (optional) initial state vector (10 x 1). Defaults to the +% operating steady state from pke_initial_conditions(plant). % % Outputs: % t - time vector (M x 1) from the ODE solver @@ -28,7 +30,9 @@ function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan) fprintf('================================\n\n'); %% ===== Solve closed-loop ===== - x0 = pke_initial_conditions(plant); + if nargin < 6 || isempty(x0) + x0 = pke_initial_conditions(plant); + end options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10); rhs = @(t, x) pke_th_rhs(t, x, plant, Q_sg, ctrl_fn(t, x, plant, ref)); diff --git a/plant-model/test_linearize.m b/plant-model/test_linearize.m new file mode 100644 index 0000000..3f99ffd --- /dev/null +++ b/plant-model/test_linearize.m @@ -0,0 +1,106 @@ +%% test_linearize.m — sanity-check the Jacobians against the nonlinear sim +% +% Simulate two parallel trajectories under a small Q_sg step (5% down): +% 1. Full nonlinear plant with u = 0 +% 2. Linear dx/dt = A*dx + B_w*dw propagated by ode45 +% Compare deviations from x_op. For a small-enough disturbance, the +% trajectories should overlay to within a few percent over a few minutes. + +clear; clc; close all; +addpath('controllers'); + +plant = pke_params(); +x_op = pke_initial_conditions(plant); + +[A, B, B_w, xs, us, Qs] = pke_linearize(plant, x_op, 0, plant.P0); + +%% ===== Eigenvalue/conditioning check ===== +eigs_A = eig(A); +fprintf('\n=== Linearization at operating point ===\n'); +fprintf(' ||A|| = %.3e\n', norm(A)); +fprintf(' ||B|| = %.3e\n', norm(B)); +fprintf(' ||B_w|| = %.3e\n', norm(B_w)); +fprintf(' stable? = %s (max Re(eig) = %.3e)\n', ... + iif(all(real(eigs_A) < 1e-8), 'YES', 'NO'), max(real(eigs_A))); +fprintf(' fastest pole ~ %.3e /s\n', max(abs(real(eigs_A)))); +fprintf(' slowest pole ~ %.3e /s\n', min(abs(real(eigs_A(real(eigs_A) < -1e-10))))); + +%% ===== Linear vs nonlinear trajectory under small Q_sg step ===== +% 5% down step at t = 30s. +dQ = -0.05 * plant.P0; +Q_sg = @(t) plant.P0 + dQ * (t >= 30); + +tspan = [0, 300]; + +% Nonlinear +opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10); +rhs_nl = @(t, x) pke_th_rhs(t, x, plant, Q_sg, 0); +[t_nl, X_nl] = ode15s(rhs_nl, tspan, x_op, opts); + +% Linear (deviation) +dQ_of_t = @(t) dQ * (t >= 30); +rhs_lin = @(t, dx) A*dx + B_w*dQ_of_t(t); +[t_lin, DX_lin] = ode15s(rhs_lin, tspan, zeros(size(x_op)), opts); + +% Reconstruct absolute linear state for plotting +X_lin = DX_lin + x_op.'; + +%% ===== Plot comparison — focus on thermal states ===== +CtoF = @(T) T*9/5 + 32; +figdir = fullfile('..', 'docs', 'figures'); +if ~exist(figdir, 'dir'), mkdir(figdir); end + +figure('Position', [100 80 1100 700], 'Name', 'Linear vs nonlinear sanity'); + +subplot(2,2,1); +plot(t_nl, X_nl(:,1), 'b-', 'LineWidth', 1.5); hold on; +plot(t_lin, X_lin(:,1), 'r--','LineWidth', 1.5); +xlabel('t [s]'); ylabel('n'); grid on; +title('Normalized power'); legend('nonlinear', 'linear', 'Location', 'best'); + +subplot(2,2,2); +plot(t_nl, CtoF(X_nl(:,8)), 'b-', 'LineWidth', 1.5); hold on; +plot(t_lin, CtoF(X_lin(:,8)), 'r--','LineWidth', 1.5); +xlabel('t [s]'); ylabel('T_f [F]'); grid on; +title('Fuel temperature'); + +subplot(2,2,3); +plot(t_nl, CtoF(X_nl(:,9)), 'b-', 'LineWidth', 1.5); hold on; +plot(t_lin, CtoF(X_lin(:,9)), 'r--','LineWidth', 1.5); +xlabel('t [s]'); ylabel('T_c [F]'); grid on; +title('Avg coolant'); + +subplot(2,2,4); +plot(t_nl, CtoF(X_nl(:,10)), 'b-', 'LineWidth', 1.5); hold on; +plot(t_lin, CtoF(X_lin(:,10)),'r--','LineWidth', 1.5); +xlabel('t [s]'); ylabel('T_{cold} [F]'); grid on; +title('Cold-leg coolant'); + +sgtitle('Linear-model sanity check: 5% Q_{sg} down-step, u = 0'); +exportgraphics(gcf, fullfile(figdir, 'linearize_sanity.png'), 'Resolution', 150); + +%% ===== Quantitative error at a few times ===== +fprintf('\n=== Max |linear - nonlinear| at sampled times ===\n'); +fprintf(' (normalized by operating-point magnitude where nonzero)\n'); +for ts = [60, 120, 300] + xi_nl = interp1(t_nl, X_nl, ts); + xi_lin = interp1(t_lin, X_lin, ts); + rel_err = abs(xi_nl - xi_lin) ./ max(abs(x_op.'), 1e-6); + fprintf(' t = %3d s: max rel err = %.2e (%s)\n', ts, max(rel_err), ... + state_name(find(rel_err == max(rel_err), 1))); +end + +%% ===== Save linearization ===== +save(fullfile('..', 'reachability', 'linearization_at_op.mat'), ... + 'A', 'B', 'B_w', 'xs', 'us', 'Qs', '-v7'); +fprintf('\nSaved A/B/B_w to ../reachability/linearization_at_op.mat\n'); + +%% ===== helpers ===== +function s = state_name(k) + names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'}; + s = names{k}; +end + +function y = iif(cond, a, b) + if cond, y = a; else, y = b; end +end