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

plant-model: add shutdown/heatup/scram controllers and LQR, linearize

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
Dane Sabo 2026-04-17 12:52:03 -04:00
parent cebf8c167a
commit d2997c2861
9 changed files with 503 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
plant-model/CLAUDE.md
View File

@ -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_<mode>(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.

57
plant-model/controllers/ctrl_heatup.m Normal file
View File

@ -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

53
plant-model/controllers/ctrl_operation_lqr.m Normal file
View File

@ -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

20
plant-model/controllers/ctrl_scram.m Normal file
View File

@ -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

21
plant-model/controllers/ctrl_shutdown.m Normal file
View File

@ -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

137
plant-model/main_mode_sweep.m Normal file
View File

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

54
plant-model/pke_linearize.m Normal file
View File

@ -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<NASGU> % 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

6
plant-model/pke_solver.m
View File

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

106
plant-model/test_linearize.m Normal file
View File

@ -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