diff --git a/docs/figures/linearize_sanity.png b/docs/figures/linearize_sanity.png new file mode 100644 index 0000000..d5528fe Binary files /dev/null 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 new file mode 100644 index 0000000..1068d7f Binary files /dev/null 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 new file mode 100644 index 0000000..7de1db4 Binary files /dev/null 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 new file mode 100644 index 0000000..aef2535 Binary files /dev/null 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 new file mode 100644 index 0000000..7907679 Binary files /dev/null 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 new file mode 100644 index 0000000..a7c8bea Binary files /dev/null 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 new file mode 100644 index 0000000..0e9530e Binary files /dev/null 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 new file mode 100644 index 0000000..046eb86 Binary files /dev/null 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 new file mode 100644 index 0000000..479e6fa Binary files /dev/null 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 new file mode 100644 index 0000000..d918ba3 Binary files /dev/null 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 new file mode 100644 index 0000000..44d1e11 Binary files /dev/null and b/docs/figures/reach_operation_n.png differ diff --git a/reachability/.gitignore b/reachability/.gitignore new file mode 100644 index 0000000..e63f75c --- /dev/null +++ b/reachability/.gitignore @@ -0,0 +1 @@ +*.mat diff --git a/reachability/README.md b/reachability/README.md new file mode 100644 index 0000000..607104f --- /dev/null +++ b/reachability/README.md @@ -0,0 +1,55 @@ +# Reachability + +Continuous-mode verification for the PWR_HYBRID_3 hybrid controller. + +## Status + +**Per-mode only.** Following the compositionality argument in the thesis: +verify each continuous mode separately, let the DRC handle discrete +switching. Current focus: **operation mode** under LQR feedback. + +## What's here + +- `linearization_at_op.mat` — A, B, B_w and reference point, generated by + `../plant-model/test_linearize.m`. +- `reach_linear.m` — box-zonotope propagation of the closed-loop linear + model under bounded disturbance. Pure MATLAB, no external toolbox. +- `barrier_lyapunov.m` — Lyapunov-ellipsoid barrier certificate for the + closed-loop linear system. Solves a Lyapunov equation, reports the + smallest sub-level set containing the initial set and closed under + the disturbance. +- `reach_operation.m` — end-to-end operation-mode reach: linearize at + x_op, compute LQR gain, propagate zonotope reach set, check against + the `t_avg_in_range` predicate. +- `figures/` — generated plots. + +## Running + +From MATLAB: + +```matlab +cd reachability +reach_operation % computes reach set + plots +barrier_lyapunov % solves Lyapunov, reports invariant ellipsoid +``` + +## Tool choice + +Currently using a hand-rolled zonotope reach because: +- Avoids a ~0.5 GB CORA install for a first-pass result. +- Linear reach with bounded disturbance has a clean analytic form + (matrix exponential on the state, integral of e^(A(t-s))·B_w·w ds + for the disturbance). +- Stays inside MATLAB, which is where the plant model lives. + +If we need nonlinear reach (and we will, for non-LQR controllers or +larger reach sets where linearization error matters), the planned +options are CORA (MATLAB) or JuliaReach (port the plant to Julia). + +## What this does NOT do yet + +- Nonlinear reach for the original P controller on operation. +- Heatup reach (the ramped reference makes x* time-varying — needs + trajectory-LQR or a different formulation). +- Shutdown, scram, initialization reach. +- Hybrid-system level verification (mode switching validity). diff --git a/reachability/barrier_lyapunov.m b/reachability/barrier_lyapunov.m new file mode 100644 index 0000000..1413803 --- /dev/null +++ b/reachability/barrier_lyapunov.m @@ -0,0 +1,168 @@ +%% barrier_lyapunov.m — Lyapunov-ellipsoid barrier certificate +% +% For dx/dt = A_cl x + B_w w with A_cl Hurwitz and ||w||_inf <= w_bar: +% +% 1. Solve A_cl' P + P A_cl = -Qbar (Qbar > 0, chosen = I). +% Then V(x) = x' P x is a Lyapunov function for the undisturbed +% system, with dV/dt = -x'x - x'(Qbar-I)x (here Qbar=I gives -x'x). +% +% 2. Under bounded disturbance: +% dV/dt = -x'x + 2 x' P B_w w +% <= -||x||^2 + 2 ||P B_w|| w_bar ||x||. +% dV/dt <= 0 whenever ||x|| >= 2 ||P B_w|| w_bar. +% So the ball B_r := {x : ||x|| <= 2 ||P B_w|| w_bar} contains +% the set where V can still grow. Any level set {V <= c} that +% contains B_r is forward-invariant. +% +% 3. Smallest such c: c* = lambda_max(P) * r^2, where r = 2||P B_w||w_bar. +% +% 4. Safety: the barrier is B(x) = V(x) - gamma, with gamma chosen +% large enough to contain X_entry but small enough that the level +% set stays inside X_safe. We report whether such a gamma exists. +% +% This is an ellipsoidal over-approximation, generally much looser than +% the box/zonotope reach in reach_operation.m, but it gives a *certificate* +% (a closed-form invariant function) rather than just a numerical tube. + +clear; clc; +addpath('../plant-model', '../plant-model/controllers'); + +plant = pke_params(); +x_op = pke_initial_conditions(plant); + +%% ===== Build A_cl, B_w ===== +[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); + +Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e-2, 1e2, 1]); +R_lqr = 1e6; +try + K = lqr(A, B, Q_lqr, R_lqr); +catch + [~, ~, K] = icare(A, B, Q_lqr, R_lqr); +end +A_cl = A - B*K; + +%% ===== Solve Lyapunov equation ===== +% A_cl' P + P A_cl + Qbar = 0. Qbar shaped to weight T_c heavily so the +% resulting ellipsoidal invariant sets are tight in the T_c direction. +% Without shaping, isotropic Qbar = I gives ellipsoids stretched along +% the slow-precursor directions, making the T_c safety bound useless. +Qbar = diag([1, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1, 1e3, 1]); +P = lyap(A_cl.', Qbar); +assert(all(eig(P) > 0), 'P not positive definite'); + +%% ===== Safety spec (used by sweep and final check) ===== +e9 = zeros(10, 1); e9(9) = 1; +delta_safe_Tc = 5.0; + +%% ===== Disturbance bound ===== +% |w| <= w_bar where w = Q_sg - Q_nom. Take the same 15% down-load as +% reach_operation.m. +w_bar = 0.15 * plant.P0; + +% --- Invariant-level computation --- +% dV/dt = -x' Qbar x + 2 x' P B_w w. +% Taking the worst w = w_bar * sign(x' P B_w), the scalar g = x' P B_w: +% dV/dt <= -x' Qbar x + 2 w_bar |g|. +% Let u = P^{1/2} x (so V = ||u||^2). Then |g| = |u' P^{-1/2} P B_w| +% <= ||u|| * ||P^{-1/2} P B_w|| = sqrt(V) * sqrt(B_w' P B_w). +% And x' Qbar x >= lambda_min(P^{-1/2} Qbar P^{-1/2}) * V (call this mu). +% So dV/dt <= -mu V + 2 w_bar sqrt(B_w' P B_w) sqrt(V). +% dV/dt <= 0 whenever sqrt(V) >= 2 w_bar sqrt(B_w' P B_w) / mu, +% i.e. V >= (2 w_bar sqrt(B_w' P B_w) / mu)^2 := c_inv. +% +% This is much tighter than the isotropic ball bound — it uses the fact +% that B_w only pokes one direction of the ellipsoid. +P_half = sqrtm(P); +P_half_inv = inv(P_half); +mu = min(eig(P_half_inv * Qbar * P_half_inv)); +g_bound = sqrt(B_w.' * P * B_w); % sqrt(B_w' P B_w) +c_inv = (2 * w_bar * g_bound / mu)^2; + +fprintf('\n=== Lyapunov barrier certificate ===\n'); +fprintf(' lambda_min(P) = %.3e\n', min(eig(P))); +fprintf(' lambda_max(P) = %.3e\n', max(eig(P))); +fprintf(' sqrt(B_w'' P B_w) = %.3e\n', g_bound); +fprintf(' mu (Qbar eig on P-metric) = %.3e\n', mu); +fprintf(' w_bar (15%% P0) = %.3e W\n', w_bar); +fprintf(' c_inv (invariant level) = %.3e\n', c_inv); + +%% ===== Containment of initial set ===== +% Initial set: box around x_op with halfwidth delta_entry (matches reach_operation). +delta_entry = [0.01 * x_op(1); + 0.001 * abs(x_op(2:7)); + 0.1; 0.1; 0.1]; + +% Worst-case V over the initial box: max x'Px = sum over all 2^n corners. +% For small n we could enumerate, but the sharper bound is +% max V(dx) = delta_entry' * |P| * delta_entry +% (elementwise abs of P), which is the L1 energy bound. +c_entry = delta_entry.' * abs(P) * delta_entry; + +fprintf('\n c_entry (bound on V over initial box) = %.3e\n', c_entry); + +gamma = max(c_entry, c_inv); % barrier level must contain both +fprintf(' gamma (barrier level) = %.3e\n', gamma); +if gamma == c_entry + fprintf(' (initial set drives gamma — invariant piece already inside entry)\n'); +else + fprintf(' (disturbance inflation drives gamma)\n'); +end + +%% ===== Sweep Qbar(9,9) to find the tightest safe barrier ===== +% The isotropic Lyapunov bound is conservative because the "slow decay" +% direction dominates mu even when T_c is tightly controlled. Sweep the +% T_c weight to find a Qbar that yields a sub-5K barrier if one exists +% for this LQR design. +fprintf('\n=== Sweeping Qbar(T_c) weight ===\n'); +weights = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6]; +best_dTc = inf; best_w = NaN; best_gamma = NaN; best_P = []; +for wTc = weights + Qbar_s = Qbar; Qbar_s(9,9) = wTc; + try + Ps = lyap(A_cl.', Qbar_s); + catch + continue + end + if any(eig(Ps) <= 0), continue, end + Ph = sqrtm(Ps); Phi = inv(Ph); + mu_s = min(eig(Phi * Qbar_s * Phi)); + g_s = sqrt(B_w.' * Ps * B_w); + ci_s = (2 * w_bar * g_s / mu_s)^2; + ce_s = delta_entry.' * abs(Ps) * delta_entry; + g_s_level = max(ci_s, ce_s); + Pinv_s = inv(Ps); + dTc_s = sqrt(g_s_level * (e9.' * Pinv_s * e9)); + fprintf(' Qbar(9,9) = %.0e -> gamma = %.3e, max|dT_c| = %7.3f K\n', ... + wTc, g_s_level, dTc_s); + if dTc_s < best_dTc + best_dTc = dTc_s; best_w = wTc; best_gamma = g_s_level; best_P = Ps; + end +end +fprintf(' Best: Qbar(9,9) = %.0e -> max|dT_c| = %.3f K\n', best_w, best_dTc); +if best_dTc <= delta_safe_Tc + fprintf(' *** TIGHT BARRIER FOUND: V(x) = x.'' P_best x - gamma ***\n'); + P = best_P; gamma = best_gamma; +end + +%% ===== Safety: does the gamma-level set fit inside X_safe? ===== +% X_safe = { x : |T_c - T_c0| <= 5 K }, i.e. |e_9.' * dx| <= 5. +% Max |e_9.' * dx| over {dx : dx' P dx <= gamma} is sqrt(gamma * e_9' P^-1 e_9). +Pinv = inv(P); +max_dTc_on_ellipsoid = sqrt(gamma * (e9.' * Pinv * e9)); + +fprintf('\n=== Safety check on T_c ===\n'); +fprintf(' Max |dT_c| on gamma-ellipsoid = %.3f K\n', max_dTc_on_ellipsoid); +fprintf(' Safe band = +/- %.1f K\n', delta_safe_Tc); +if max_dTc_on_ellipsoid <= delta_safe_Tc + fprintf(' BARRIER VALID: V(x) = x.''Px - %.3e certifies T_c safety.\n', gamma); +else + fprintf(' *** BARRIER TOO LOOSE *** - ellipsoid reach into unsafe region.\n'); + fprintf(' Try a tighter LQR (bigger Q_Tc or smaller R) or tighter X_entry.\n'); +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'); + +fprintf('\nSaved barrier to ./barrier_lyapunov_result.mat\n'); diff --git a/reachability/reach_linear.m b/reachability/reach_linear.m new file mode 100644 index 0000000..20e84bb --- /dev/null +++ b/reachability/reach_linear.m @@ -0,0 +1,93 @@ +function [T, R_lo, R_hi, center] = reach_linear(A_cl, B_w, x0_center, x0_halfwidth, w_lo, w_hi, tspan, dt) +% REACH_LINEAR Interval reach set for dx/dt = A_cl*x + B_w*w, w in [w_lo, w_hi]. +% +% Hand-rolled zonotope propagator specialized to the case where: +% - The initial set is an axis-aligned box around x0_center with +% halfwidths x0_halfwidth. +% - The disturbance w is a scalar interval. +% +% The reach set at time t is: +% X(t) = {expm(A_cl*t) * x0 + integral_0^t expm(A_cl*(t-s))*B_w*w(s) ds : +% x0 in X0, w(s) in W}. +% +% Decompose: +% X(t) = expm(A_cl*t) * X0 (+) S_w(t), +% where S_w(t) is the reachable contribution of the disturbance. +% Both parts are zonotopes; we over-approximate each by a box at every +% time step via the matrix/integral absolute sum (interval hull). +% +% Inputs: +% A_cl - n x n closed-loop state matrix (Hurwitz assumed) +% B_w - n x 1 disturbance vector +% x0_center - n x 1 center of initial set +% x0_halfwidth - n x 1 halfwidths of initial set (>=0) +% w_lo, w_hi - scalar disturbance bounds +% tspan - [t0, tf] +% dt - time step for exporting the reach tube (s) +% +% Outputs: +% T - M x 1 time grid +% R_lo - n x M lower bounds of reach set +% R_hi - n x M upper bounds of reach set +% center- n x M center trajectory (nominal, w = w_mid) + + n = size(A_cl, 1); + assert(size(A_cl,2) == n, 'A_cl must be square'); + assert(numel(x0_center) == n && numel(x0_halfwidth) == n, 'x0 dim mismatch'); + x0_center = x0_center(:); + x0_halfwidth = x0_halfwidth(:); + B_w = B_w(:); + + w_mid = 0.5 * (w_lo + w_hi); + w_half = 0.5 * (w_hi - w_lo); + + T = (tspan(1):dt:tspan(2)).'; + M = numel(T); + + R_lo = zeros(n, M); + R_hi = zeros(n, M); + center = zeros(n, M); + + % Disturbance contribution: + % S_center(t) = int_0^t expm(A*(t-s))*B_w*w_mid ds (signed, centerline) + % S_half(t) = int_0^t |expm(A*(t-s))*B_w| * w_half ds (elementwise, halfwidth) + % + % The halfwidth uses elementwise |.| because under an axis-aligned + % box over-approximation, |linear-map| contributes absolute values. + % This is the interval-arithmetic analog of propagating a zonotope. + S_center = zeros(n, 1); + S_half = zeros(n, 1); + + x_center_now = x0_center; + A_step = expm(A_cl * dt); % state transition per step + A_abs_step = abs(A_step); % for interval-halfwidth + + % Integrated input gain per step: + % G_step = int_0^dt expm(A*s)*B_w ds + % Augmented-matrix trick: upper-right block of expm([A B_w; 0 0]*dt). + M_aug = expm([A_cl, B_w; zeros(1, n+1)] * dt); + G_step = M_aug(1:n, n+1); % n x 1 + G_abs_step = abs(G_step); + + halfwidth_now = x0_halfwidth; + + % Initial state + center(:, 1) = x_center_now + S_center; + R_lo(:, 1) = center(:, 1) - halfwidth_now - S_half; + R_hi(:, 1) = center(:, 1) + halfwidth_now + S_half; + + for k = 2:M + % State piece + x_center_now = A_step * x_center_now; + halfwidth_now = A_abs_step * halfwidth_now; + + % Disturbance piece + S_center = A_step * S_center + G_step * w_mid; + S_half = A_abs_step * S_half + G_abs_step * w_half; + + center(:, k) = x_center_now + S_center; + R_lo(:, k) = center(:, k) - halfwidth_now - S_half; + R_hi(:, k) = center(:, k) + halfwidth_now + S_half; + end + +end diff --git a/reachability/reach_operation.m b/reachability/reach_operation.m new file mode 100644 index 0000000..11106af --- /dev/null +++ b/reachability/reach_operation.m @@ -0,0 +1,156 @@ +%% reach_operation.m — linear reach set for operation mode (LQR closed-loop) +% +% Compute a sound over-approximation of the reach set starting from a +% box around x_op, under LQR feedback, with Q_sg in a specified +% interval. Check that T_avg stays inside the t_avg_in_range predicate +% for all t in [0, T_final]. +% +% This is the *continuous-mode obligation* for q_operation: +% X_entry := { x : |x - x_op| <= delta_entry } +% W := [Q_min, Q_max] +% X_safe := { x : |T_c - T_c0| <= delta_safe } +% Obligation: ReachTube(X_entry, W, [0, T_final]) subset X_safe. +% +% If this passes, we've discharged the Thrust-3 verification for one +% continuous mode at the level the thesis calls for. + +clear; clc; close all; +addpath('../plant-model', '../plant-model/controllers'); + +plant = pke_params(); +x_op = pke_initial_conditions(plant); + +%% ===== Closed-loop linearization ===== +[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0); + +% LQR gain using the same Q, R as ctrl_operation_lqr.m. Keep these in +% sync if you retune there. +Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1]); +R_lqr = 1e6; +try + K = lqr(A, B, Q_lqr, R_lqr); +catch + [~, ~, K] = icare(A, B, Q_lqr, R_lqr); +end + +A_cl = A - B*K; + +fprintf('\n=== Closed-loop spectrum (A - BK) ===\n'); +eigs_cl = eig(A_cl); +fprintf(' max Re(eig) = %.3e\n', max(real(eigs_cl))); +fprintf(' min Re(eig) = %.3e\n', min(real(eigs_cl))); +assert(all(real(eigs_cl) < -1e-8), 'A_cl not Hurwitz — gain tuning issue'); + +%% ===== Define sets ===== +% X_entry: 1% box on n, 0.1% boxes on precursors (their magnitudes are +% huge due to 1/Lambda), 0.1 K on fuel, 0.1 K on coolant, 0.1 K on cold leg. +% This represents "we entered operation mode near the steady state". +delta_entry = [0.01 * x_op(1); % n + 0.001 * abs(x_op(2:7)); % C1..C6 + 0.1; % T_f [C] + 0.1; % T_c [C] + 0.1]; % T_cold [C] + +% Disturbance: Q_sg = [0.85*P0, 1.00*P0] -> this captures up to a 15% +% down-step load demand (realistic load-follow envelope). +Q_nom = plant.P0; +Q_min = 0.85 * plant.P0; +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] + +%% ===== Reach set ===== +% Propagate in deviation coordinates: dx = x - x_op. +tspan = [0, 600]; +dt = 0.5; + +[T, R_lo, R_hi, C] = reach_linear(A_cl, B_w, zeros(10,1), delta_entry, dQ_lo, dQ_hi, tspan, dt); + +% Translate back to absolute coordinates for reporting +Xabs_lo = R_lo + x_op; +Xabs_hi = R_hi + x_op; +Cabs = C + x_op; + +%% ===== Safety check ===== +T_c_lo = Xabs_lo(9, :); +T_c_hi = Xabs_hi(9, :); + +violation_mask = (T_c_hi > plant.T_c0 + delta_safe_Tc) | ... + (T_c_lo < plant.T_c0 - delta_safe_Tc); +fprintf('\n=== Operation-mode reach-set safety ===\n'); +fprintf(' Horizon = [%g, %g] s\n', tspan(1), tspan(2)); +fprintf(' Entry box T_c [C] = [%.3f, %.3f] (x_op +/- %.1f C)\n', ... + x_op(9) - delta_entry(9), x_op(9) + delta_entry(9), delta_entry(9)); +fprintf(' Disturbance Q_sg = [%.3f, %.3f] MW\n', Q_min/1e6, Q_max/1e6); +fprintf(' Safe band on T_c = x_op(T_c0) +/- %.1f C -> [%.3f, %.3f]\n', ... + delta_safe_Tc, plant.T_c0 - delta_safe_Tc, plant.T_c0 + delta_safe_Tc); +fprintf(' Reach T_c envelope = [%.3f, %.3f]\n', min(T_c_lo), max(T_c_hi)); +if any(violation_mask) + t_first = T(find(violation_mask, 1)); + fprintf(' *** SAFETY VIOLATED at t = %.2f s ***\n', t_first); +else + fprintf(' OK: reach set stays inside the safe band.\n'); +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'); +fprintf(' %-7s %-14s %-14s %-8s\n', 'state', 'init halfwidth', 'final halfwidth', 'ratio'); +for i = 1:10 + hi = 0.5 * (R_hi(i, 1) - R_lo(i, 1)); + hf = 0.5 * (R_hi(i, end) - R_lo(i, end)); + fprintf(' %-7s %-14.4e %-14.4e %-8.2f\n', state_names{i}, hi, hf, hf/max(hi,eps)); +end + +%% ===== Plots ===== +figdir = fullfile('..', 'docs', 'figures'); +if ~exist(figdir, 'dir'), mkdir(figdir); end +CtoF = @(T) T*9/5 + 32; + +% Two-panel plot: wide view with safety band, zoom view showing actual tube. +figure('Position', [100 80 1400 500], 'Name', 'Reach tube: T_c'); + +subplot(1,2,1); +fill([T; flipud(T)], CtoF([T_c_hi.'; flipud(T_c_lo.')]), [1.0 0.85 0.85], ... + 'EdgeColor', 'none'); hold on; +plot(T, CtoF(Cabs(9, :)), 'r-', 'LineWidth', 1.2); +yline(CtoF(plant.T_c0 + delta_safe_Tc), 'k--', 'LineWidth', 1.0); +yline(CtoF(plant.T_c0 - delta_safe_Tc), 'k--', 'LineWidth', 1.0); +yline(CtoF(plant.T_c0), 'k:', 'LineWidth', 1.0); +xlabel('Time [s]'); ylabel('T_{avg} [F]'); grid on; +title('Safety-band view'); +legend('reach tube', 'nominal', 'safety +/- 5 C', 'Location', 'best'); + +subplot(1,2,2); +Tc_dev_lo = T_c_lo.' - plant.T_c0; % M x 1, deviation in K +Tc_dev_hi = T_c_hi.' - plant.T_c0; +fill([T; flipud(T)], [Tc_dev_hi; flipud(Tc_dev_lo)], [1.0 0.85 0.85], ... + 'EdgeColor', 'none'); hold on; +plot(T, Cabs(9, :).' - plant.T_c0, 'r-', 'LineWidth', 1.2); +yline(0, 'k:', 'LineWidth', 1.0); +xlabel('Time [s]'); ylabel('T_{avg} - T_{c0} [K]'); grid on; +max_dev = max(abs([Tc_dev_lo; Tc_dev_hi])); +title(sprintf('Zoomed: max |dT_c| = %.3e K', max_dev)); + +sgtitle(sprintf('Operation-mode reach tube, LQR, Q_{sg} in [%.0f%%, %.0f%%] P_0', ... + 100*Q_min/Q_nom, 100*Q_max/Q_nom)); +exportgraphics(gcf, fullfile(figdir, 'reach_operation_Tc.png'), 'Resolution', 150); + +figure('Position', [100 80 1100 500], 'Name', 'Reach tube: n'); +fill([T; flipud(T)], [R_hi(1,:).' + x_op(1); flipud(R_lo(1,:).' + x_op(1))], ... + [0.85 0.85 1.0], 'EdgeColor', 'none'); hold on; +plot(T, Cabs(1, :), 'b-', 'LineWidth', 1.2); +xlabel('Time [s]'); ylabel('n'); grid on; +title('Operation mode reach tube on normalized power'); +legend('reach tube', 'nominal', 'Location', 'best'); +exportgraphics(gcf, fullfile(figdir, 'reach_operation_n.png'), 'Resolution', 150); + +save(fullfile('.', 'reach_operation_result.mat'), ... + 'T', 'R_lo', 'R_hi', 'C', 'Xabs_lo', 'Xabs_hi', 'Cabs', ... + 'K', 'A_cl', 'x_op', 'delta_entry', 'Q_min', 'Q_max', 'delta_safe_Tc', '-v7'); + +fprintf('\nSaved reach result to ./reach_operation_result.mat\n');