diff --git a/code/scripts/plot_reach_tubes.jl b/code/scripts/plot_reach_tubes.jl new file mode 100644 index 0000000..76de00c --- /dev/null +++ b/code/scripts/plot_reach_tubes.jl @@ -0,0 +1,168 @@ +#!/usr/bin/env julia +# +# plot_reach_tubes.jl — multi-panel reach-tube visualization. +# +# Produces a four-panel figure from a reach-result .mat file: +# (1) Temperature tube overlay: T_c, T_hot, T_cold envelopes together. +# The gap between T_hot and T_cold is a direct proxy for core +# thermal power (P = W·c_c·ΔT). +# (2) ΔT_core = T_hot - T_cold envelope (the power proxy, explicit). +# (3) Reactivity ρ envelope, normalized by β (in dollars). +# (4) Normalized power n envelope. +# +# Two input formats supported: +# operation: reach_operation_result.mat (linear reach, has R_lo/R_hi +# matrices, time vector T, nominal X_nom). +# heatup_pj: reach_heatup_pj_tight_full.mat (per-timestep envelopes +# Tc_lo_ts/Tc_hi_ts/... already extracted; rho from PJ +# controller). +# +# Usage: +# julia --project=. scripts/plot_reach_tubes.jl operation +# julia --project=. scripts/plot_reach_tubes.jl heatup_pj + +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) + +using MAT +using Plots +gr() + +include(joinpath(@__DIR__, "..", "src", "pke_params.jl")) +const PLANT = pke_params() + +function plot_tubes_operation() + mat_path = joinpath(@__DIR__, "..", "..", "reachability", "reach_operation_result.mat") + d = matread(mat_path) + + T_vec = vec(d["T"]) # time grid + X_lo = d["X_lo"] # 10 × M lower bounds + X_hi = d["X_hi"] # 10 × M upper bounds + X_nom = d["X_nom"] # 10 × M nominal + + # States at indices: n=1, T_f=8, T_c=9, T_cold=10. + n_lo = X_lo[1, :]; n_hi = X_hi[1, :] + Tf_lo = X_lo[8, :]; Tf_hi = X_hi[8, :] + Tc_lo = X_lo[9, :]; Tc_hi = X_hi[9, :] + Tco_lo = X_lo[10, :]; Tco_hi = X_hi[10, :] + + # T_hot = 2*T_c - T_cold; envelope via worst-case signed combination. + Th_lo = 2 .* Tc_lo .- Tco_hi + Th_hi = 2 .* Tc_hi .- Tco_lo + + # ΔT_core = T_hot - T_cold = 2*(T_c - T_cold). + dT_lo = 2 .* (Tc_lo .- Tco_hi) + dT_hi = 2 .* (Tc_hi .- Tco_lo) + + # ρ under LQR: ρ_total = u + α_f·dT_f + α_c·dT_c where u = -K(x-x_op). + # For a quick envelope, compute ρ at the nominal trajectory and + # inflate by bounds of the feedback contributions from the tube. + # Simpler: use total reactivity = u + α_f*(T_f-T_f0) + α_c*(T_c-T_c0). + # u under LQR is small; we approximate ρ envelope by α feedback + # alone (the u contribution is ≤ few cents). + rho_nom = PLANT.alpha_f .* (X_nom[8, :] .- PLANT.T_f0) .+ + PLANT.alpha_c .* (X_nom[9, :] .- PLANT.T_c0) + # Envelope via worst-case T_f, T_c. + rho_lo = PLANT.alpha_f .* (Tf_hi .- PLANT.T_f0) .+ + PLANT.alpha_c .* (Tc_hi .- PLANT.T_c0) # both α < 0, max T → min ρ + rho_hi = PLANT.alpha_f .* (Tf_lo .- PLANT.T_f0) .+ + PLANT.alpha_c .* (Tc_lo .- PLANT.T_c0) + + title_stem = "Operation-mode LQR reach tubes" + _plot_common(T_vec, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi, + dT_lo, dT_hi, rho_lo, rho_hi, n_lo, n_hi, title_stem, + "reach_operation_tubes.png") +end + +function plot_tubes_heatup_pj() + mat_path = joinpath(@__DIR__, "..", "..", "reachability", + "reach_heatup_pj_tight_full.mat") + d = matread(mat_path) + + t_arr = vec(d["t_arr"]) + Tc_lo = vec(d["Tc_lo_ts"]); Tc_hi = vec(d["Tc_hi_ts"]) + Tf_lo = vec(d["Tf_lo_ts"]); Tf_hi = vec(d["Tf_hi_ts"]) + Tco_lo = vec(d["Tco_lo_ts"]); Tco_hi = vec(d["Tco_hi_ts"]) + n_lo = vec(d["n_lo_ts"]); n_hi = vec(d["n_hi_ts"]) + rho_lo = vec(d["rho_lo_ts"]); rho_hi = vec(d["rho_hi_ts"]) + + Th_lo = 2 .* Tc_lo .- Tco_hi + Th_hi = 2 .* Tc_hi .- Tco_lo + dT_lo = 2 .* (Tc_lo .- Tco_hi) + dT_hi = 2 .* (Tc_hi .- Tco_lo) + + title_stem = "Heatup PJ (tight entry) reach tubes" + _plot_common(t_arr, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi, + dT_lo, dT_hi, rho_lo, rho_hi, n_lo, n_hi, title_stem, + "reach_heatup_pj_tubes.png") +end + +function _plot_common(t, Tc_lo, Tc_hi, Th_lo, Th_hi, Tco_lo, Tco_hi, + dT_lo, dT_hi, rho_lo, rho_hi, n_lo, n_hi, + title_stem, outname) + CtoF(T) = T*9/5 + 32 + + # Panel 1: T_c / T_hot / T_cold overlaid. + p1 = plot(xlabel="Time [s]", ylabel="T [°C]", + title="T tubes", legend=:right) + plot!(p1, t, Tc_hi, fillrange=Tc_lo, fillalpha=0.30, color=:red, + linealpha=0, label="T_c tube") + plot!(p1, t, Th_hi, fillrange=Th_lo, fillalpha=0.25, color=:orange, + linealpha=0, label="T_hot tube") + plot!(p1, t, Tco_hi, fillrange=Tco_lo, fillalpha=0.25, color=:blue, + linealpha=0, label="T_cold tube") + + # Panel 2: ΔT_core = T_hot - T_cold (power proxy at constant flow). + P_lo_MW = PLANT.W * PLANT.c_c .* dT_lo ./ 1e6 + P_hi_MW = PLANT.W * PLANT.c_c .* dT_hi ./ 1e6 + p2 = plot(xlabel="Time [s]", ylabel="ΔT_core = T_hot - T_cold [K]", + title="Core ΔT envelope (power proxy)", legend=:right) + plot!(p2, t, dT_hi, fillrange=dT_lo, fillalpha=0.30, color=:purple, + linealpha=0, label="ΔT_core [K]") + p2b = twinx(p2) + plot!(p2b, t, P_hi_MW, fillrange=P_lo_MW, fillalpha=0.0, color=:gray, + linealpha=0.5, linestyle=:dot, label="P=W·c_c·ΔT [MW]", + ylabel="Primary thermal power [MWth]") + + # Panel 3: ρ envelope in dollars. + rho_lo_d = rho_lo ./ PLANT.beta + rho_hi_d = rho_hi ./ PLANT.beta + p3 = plot(xlabel="Time [s]", ylabel="ρ [\$]", + title="Reactivity envelope (1 \$ = β = prompt-critical)", + legend=:right) + plot!(p3, t, rho_hi_d, fillrange=rho_lo_d, fillalpha=0.3, + color=:darkgreen, linealpha=0, label="ρ / β") + hline!(p3, [1.0, -1.0], ls=:dash, color=:red, + label="prompt ±1 \$") + hline!(p3, [0.0], ls=:dot, color=:black, label="critical") + + # Panel 4: n envelope (log scale if needed). + p4 = plot(xlabel="Time [s]", ylabel="n (normalized power)", + title="n envelope", legend=:right) + plot!(p4, t, n_hi, fillrange=n_lo, fillalpha=0.3, color=:black, + linealpha=0, label="n tube") + + fig = plot(p1, p2, p3, p4, layout=(2, 2), size=(1300, 800), + plot_title=title_stem) + figdir = joinpath(@__DIR__, "..", "..", "docs", "figures") + isdir(figdir) || mkpath(figdir) + outpath = joinpath(figdir, outname) + savefig(fig, outpath) + println("Saved $outpath") +end + +# CLI dispatch. +which_plot = length(ARGS) > 0 ? ARGS[1] : "both" +if which_plot in ("operation", "both") + plot_tubes_operation() +end +if which_plot in ("heatup_pj", "both") + mat_path = joinpath(@__DIR__, "..", "..", "reachability", + "reach_heatup_pj_tight_full.mat") + if isfile(mat_path) + plot_tubes_heatup_pj() + else + println("Skipping heatup_pj plot — $mat_path not found.") + println("(Run scripts/reach_heatup_pj_tight_full.jl first.)") + end +end diff --git a/code/scripts/reach_heatup_pj.jl b/code/scripts/reach_heatup_pj.jl index bf49613..5ddd4bb 100644 --- a/code/scripts/reach_heatup_pj.jl +++ b/code/scripts/reach_heatup_pj.jl @@ -138,35 +138,41 @@ for T_probe in probe_horizons println(" TMJets: $n_sets reach-sets, wall-time $(round(elapsed; digits=1)) s") flow_hr = overapproximate(flow, Hyperrectangle) - # Envelope at the FINAL set. - Tc_lo_env = minimum(low(set(R), 8) for R in flow_hr) - Tc_hi_env = maximum(high(set(R), 8) for R in flow_hr) - Tf_lo_env = minimum(low(set(R), 7) for R in flow_hr) - Tf_hi_env = maximum(high(set(R), 7) for R in flow_hr) - Tcold_lo = minimum(low(set(R), 9) for R in flow_hr) - Tcold_hi = maximum(high(set(R), 9) for R in flow_hr) - # Reconstruct n envelope at each step from C and T_c via PJ formula. - n_env_lo = Inf - n_env_hi = -Inf - for R in flow_hr + n_steps = length(flow_hr) + # Per-timestep envelopes for plotting (saved below). + t_arr = zeros(n_steps) + Tc_lo_ts = zeros(n_steps); Tc_hi_ts = zeros(n_steps) + Tf_lo_ts = zeros(n_steps); Tf_hi_ts = zeros(n_steps) + Tco_lo_ts = zeros(n_steps); Tco_hi_ts = zeros(n_steps) + n_lo_ts = zeros(n_steps); n_hi_ts = zeros(n_steps) + rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps) + for (k, R) in enumerate(flow_hr) s = set(R) + t_arr[k] = high(s, 10) # time-state upper bound ≈ current time + Tc_lo_ts[k] = low(s, 8); Tc_hi_ts[k] = high(s, 8) + Tf_lo_ts[k] = low(s, 7); Tf_hi_ts[k] = high(s, 7) + Tco_lo_ts[k] = low(s, 9); Tco_hi_ts[k] = high(s, 9) + # Algebraic n reconstruction at this reach-set. sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) + LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6) sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) + LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6) - # rho range: ρ = Kp*(T_ref - T_c). T_ref bounded by [T_STANDBY, T_C0], - # T_c bounded by current envelope. - rho_lo = KP_HEATUP * (T_STANDBY - high(s, 8)) - rho_hi = KP_HEATUP * (T_C0 - low(s, 8)) - denom_lo = BETA - rho_hi # smaller denom => larger n - denom_hi = BETA - rho_lo + rho_lo_here = KP_HEATUP * (T_STANDBY - high(s, 8)) + rho_hi_here = KP_HEATUP * (T_C0 - low(s, 8)) + rho_lo_ts[k] = rho_lo_here + rho_hi_ts[k] = rho_hi_here + denom_lo = BETA - rho_hi_here + denom_hi = BETA - rho_lo_here if denom_lo > 0 - n_lo_local = LAMBDA * sumLC_lo / denom_hi - n_hi_local = LAMBDA * sumLC_hi / denom_lo - n_env_lo = min(n_env_lo, n_lo_local) - n_env_hi = max(n_env_hi, n_hi_local) + n_lo_ts[k] = LAMBDA * sumLC_lo / denom_hi + n_hi_ts[k] = LAMBDA * sumLC_hi / denom_lo end end + # Also global-envelope values (for backward compat). + Tc_lo_env = minimum(Tc_lo_ts); Tc_hi_env = maximum(Tc_hi_ts) + Tf_lo_env = minimum(Tf_lo_ts); Tf_hi_env = maximum(Tf_hi_ts) + Tcold_lo = minimum(Tco_lo_ts); Tcold_hi = maximum(Tco_hi_ts) + n_env_lo = minimum(n_lo_ts); n_env_hi = maximum(n_hi_ts) println(" n envelope (reconstructed): [$(round(n_env_lo; sigdigits=4)), $(round(n_env_hi; sigdigits=4))]") println(" T_f envelope: [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C") @@ -176,7 +182,13 @@ for T_probe in probe_horizons Tc=(Tc_lo_env, Tc_hi_env), Tf=(Tf_lo_env, Tf_hi_env), Tcold=(Tcold_lo, Tcold_hi), - n=(n_env_lo, n_env_hi)) + n=(n_env_lo, n_env_hi), + t_arr=t_arr, + Tc_lo_ts=Tc_lo_ts, Tc_hi_ts=Tc_hi_ts, + Tf_lo_ts=Tf_lo_ts, Tf_hi_ts=Tf_hi_ts, + Tco_lo_ts=Tco_lo_ts, Tco_hi_ts=Tco_hi_ts, + n_lo_ts=n_lo_ts, n_hi_ts=n_hi_ts, + rho_lo_ts=rho_lo_ts, rho_hi_ts=rho_hi_ts) catch err msg = sprint(showerror, err) println(" FAILED: ", first(msg, 300)) @@ -204,14 +216,27 @@ for T_probe in probe_horizons haskey(results, T_probe) || continue r = results[T_probe] if r.status == "OK" - saved["T_$(Int(T_probe))_Tc_lo"] = r.Tc[1] - saved["T_$(Int(T_probe))_Tc_hi"] = r.Tc[2] - saved["T_$(Int(T_probe))_Tf_lo"] = r.Tf[1] - saved["T_$(Int(T_probe))_Tf_hi"] = r.Tf[2] - saved["T_$(Int(T_probe))_Tcold_lo"] = r.Tcold[1] - saved["T_$(Int(T_probe))_Tcold_hi"] = r.Tcold[2] - saved["T_$(Int(T_probe))_n_lo"] = r.n[1] - saved["T_$(Int(T_probe))_n_hi"] = r.n[2] + pre = "T_$(Int(T_probe))_" + saved[pre * "Tc_lo"] = r.Tc[1] + saved[pre * "Tc_hi"] = r.Tc[2] + saved[pre * "Tf_lo"] = r.Tf[1] + saved[pre * "Tf_hi"] = r.Tf[2] + saved[pre * "Tcold_lo"] = r.Tcold[1] + saved[pre * "Tcold_hi"] = r.Tcold[2] + saved[pre * "n_lo"] = r.n[1] + saved[pre * "n_hi"] = r.n[2] + # Per-timestep envelopes + saved[pre * "t_arr"] = r.t_arr + saved[pre * "Tc_lo_ts"] = r.Tc_lo_ts + saved[pre * "Tc_hi_ts"] = r.Tc_hi_ts + saved[pre * "Tf_lo_ts"] = r.Tf_lo_ts + saved[pre * "Tf_hi_ts"] = r.Tf_hi_ts + saved[pre * "Tco_lo_ts"] = r.Tco_lo_ts + saved[pre * "Tco_hi_ts"] = r.Tco_hi_ts + saved[pre * "n_lo_ts"] = r.n_lo_ts + saved[pre * "n_hi_ts"] = r.n_hi_ts + saved[pre * "rho_lo_ts"] = r.rho_lo_ts + saved[pre * "rho_hi_ts"] = r.rho_hi_ts end end matwrite(mat_out, saved) diff --git a/code/scripts/reach_heatup_pj_tight_full.jl b/code/scripts/reach_heatup_pj_tight_full.jl new file mode 100644 index 0000000..5808efe --- /dev/null +++ b/code/scripts/reach_heatup_pj_tight_full.jl @@ -0,0 +1,145 @@ +#!/usr/bin/env julia +# +# reach_heatup_pj_tight_full.jl — tight-entry heatup PJ reach, +# per-timestep envelopes saved for plotting (T_c, T_h, T_cold, rho, n). +# +# Identical dynamics to reach_heatup_pj_tight.jl but keeps the full +# tube data so we can overlay T_c and T_h (and infer power) on one plot. + +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) + +using LinearAlgebra +using ReachabilityAnalysis, LazySets +using MAT + +const LAMBDA = 1e-4 +const BETA_1, BETA_2, BETA_3, BETA_4, BETA_5, BETA_6 = + 0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273 +const BETA = BETA_1 + BETA_2 + BETA_3 + BETA_4 + BETA_5 + BETA_6 +const LAM_1, LAM_2, LAM_3, LAM_4, LAM_5, LAM_6 = + 0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01 +const P0 = 1e9 +const M_F, C_F, M_C, C_C, HA, W_M, M_SG = + 50000.0, 300.0, 20000.0, 5450.0, 5e7, 5000.0, 30000.0 +const T_COLD0 = 290.0 +const DT_CORE = P0 / (W_M * C_C) +const T_HOT0 = T_COLD0 + DT_CORE +const T_C0 = (T_HOT0 + T_COLD0) / 2 +const T_F0 = T_C0 + P0 / HA + +const T_STANDBY = T_C0 - 33.333333 +const RAMP_RATE_CS = 28.0 / 3600 +const KP_HEATUP = 1e-4 + +@taylorize function rhs_heatup_pj_tf!(dx, x, p, t) + rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[10] - x[8]) + sum_lam_C = LAM_1*x[1] + LAM_2*x[2] + LAM_3*x[3] + + LAM_4*x[4] + LAM_5*x[5] + LAM_6*x[6] + denom = BETA - rho + n = LAMBDA * sum_lam_C / denom + inv_factor = sum_lam_C / denom + dx[1] = BETA_1 * inv_factor - LAM_1 * x[1] + dx[2] = BETA_2 * inv_factor - LAM_2 * x[2] + dx[3] = BETA_3 * inv_factor - LAM_3 * x[3] + dx[4] = BETA_4 * inv_factor - LAM_4 * x[4] + dx[5] = BETA_5 * inv_factor - LAM_5 * x[5] + dx[6] = BETA_6 * inv_factor - LAM_6 * x[6] + dx[7] = (P0 * n - HA * (x[7] - x[8])) / (M_F * C_F) + dx[8] = (HA * (x[7] - x[8]) - 2 * W_M * C_C * (x[8] - x[9])) / (M_C * C_C) + dx[9] = (2 * W_M * C_C * (x[8] - x[9])) / (M_SG * C_C) + dx[10] = one(x[1]) + return nothing +end + +# Tight X_entry. +n_lo, n_hi = 1.0e-3, 2.0e-3 +T_f_lo, T_f_hi = 285.0, 291.0 +T_c_lo, T_c_hi = 285.0, 291.0 +T_cold_lo, T_cold_hi = 278.0, 285.0 + +n_mid = 0.5 * (n_lo + n_hi) +C_mid = [BETA_1/(LAM_1*LAMBDA), BETA_2/(LAM_2*LAMBDA), + BETA_3/(LAM_3*LAMBDA), BETA_4/(LAM_4*LAMBDA), + BETA_5/(LAM_5*LAMBDA), BETA_6/(LAM_6*LAMBDA)] .* n_mid + +x_lo = [C_mid[1]*(n_lo/n_mid), C_mid[2]*(n_lo/n_mid), + C_mid[3]*(n_lo/n_mid), C_mid[4]*(n_lo/n_mid), + C_mid[5]*(n_lo/n_mid), C_mid[6]*(n_lo/n_mid), + T_f_lo, T_c_lo, T_cold_lo, 0.0] +x_hi = [C_mid[1]*(n_hi/n_mid), C_mid[2]*(n_hi/n_mid), + C_mid[3]*(n_hi/n_mid), C_mid[4]*(n_hi/n_mid), + C_mid[5]*(n_hi/n_mid), C_mid[6]*(n_hi/n_mid), + T_f_hi, T_c_hi, T_cold_hi, 0.0] +X0 = Hyperrectangle(low=x_lo, high=x_hi) + +println("\n=== Tight-entry heatup PJ reach — saving full per-step envelopes ===") + +T_probe = 300.0 +sys = BlackBoxContinuousSystem(rhs_heatup_pj_tf!, 10) +prob = InitialValueProblem(sys, X0) +alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000) +t_start = time() +sol = solve(prob; T=T_probe, alg=alg) +println(" Wall time: $(round(time() - t_start; digits=1)) s") +flow_hr = overapproximate(flowpipe(sol), Hyperrectangle) +n_steps = length(flow_hr) +println(" $n_steps reach-sets") + +t_arr = zeros(n_steps) +Tc_lo_ts = zeros(n_steps); Tc_hi_ts = zeros(n_steps) +Tf_lo_ts = zeros(n_steps); Tf_hi_ts = zeros(n_steps) +Tco_lo_ts = zeros(n_steps); Tco_hi_ts = zeros(n_steps) +n_lo_ts = zeros(n_steps); n_hi_ts = zeros(n_steps) +rho_lo_ts = zeros(n_steps); rho_hi_ts = zeros(n_steps) + +for (k, R) in enumerate(flow_hr) + s = set(R) + t_arr[k] = high(s, 10) + Tc_lo_ts[k] = low(s, 8); Tc_hi_ts[k] = high(s, 8) + Tf_lo_ts[k] = low(s, 7); Tf_hi_ts[k] = high(s, 7) + Tco_lo_ts[k] = low(s, 9); Tco_hi_ts[k] = high(s, 9) + sumLC_lo = LAM_1*low(s,1) + LAM_2*low(s,2) + LAM_3*low(s,3) + + LAM_4*low(s,4) + LAM_5*low(s,5) + LAM_6*low(s,6) + sumLC_hi = LAM_1*high(s,1) + LAM_2*high(s,2) + LAM_3*high(s,3) + + LAM_4*high(s,4) + LAM_5*high(s,5) + LAM_6*high(s,6) + # Ramp-ref bounds at this reach-set's t. + t_hi_here = high(s, 10) + t_lo_here = low(s, 10) + # T_ref at these times (monotone increasing): + Tref_lo = T_STANDBY + RAMP_RATE_CS * t_lo_here + Tref_hi = T_STANDBY + RAMP_RATE_CS * t_hi_here + rho_lo_here = KP_HEATUP * (Tref_lo - high(s, 8)) + rho_hi_here = KP_HEATUP * (Tref_hi - low(s, 8)) + rho_lo_ts[k] = rho_lo_here + rho_hi_ts[k] = rho_hi_here + denom_lo = BETA - rho_hi_here + denom_hi = BETA - rho_lo_here + if denom_lo > 0 + n_lo_ts[k] = LAMBDA * sumLC_lo / denom_hi + n_hi_ts[k] = LAMBDA * sumLC_hi / denom_lo + end +end + +println(" T_c envelope: [$(round(minimum(Tc_lo_ts); digits=2)), $(round(maximum(Tc_hi_ts); digits=2))] °C") +println(" T_hot envelope: [$(round(minimum(2 .* Tc_lo_ts .- Tco_hi_ts); digits=2)), $(round(maximum(2 .* Tc_hi_ts .- Tco_lo_ts); digits=2))] °C") +println(" T_cold env: [$(round(minimum(Tco_lo_ts); digits=2)), $(round(maximum(Tco_hi_ts); digits=2))] °C") +println(" rho envelope: [$(round(minimum(rho_lo_ts); sigdigits=4)), $(round(maximum(rho_hi_ts); sigdigits=4))]") +println(" rho / beta: [$(round(minimum(rho_lo_ts)/BETA; digits=3)), $(round(maximum(rho_hi_ts)/BETA; digits=3))]") + +mat_out = joinpath(@__DIR__, "..", "..", "reachability", "reach_heatup_pj_tight_full.mat") +matwrite(mat_out, Dict( + "t_arr" => t_arr, + "Tc_lo_ts" => Tc_lo_ts, "Tc_hi_ts" => Tc_hi_ts, + "Tf_lo_ts" => Tf_lo_ts, "Tf_hi_ts" => Tf_hi_ts, + "Tco_lo_ts" => Tco_lo_ts, "Tco_hi_ts" => Tco_hi_ts, + "n_lo_ts" => n_lo_ts, "n_hi_ts" => n_hi_ts, + "rho_lo_ts" => rho_lo_ts, "rho_hi_ts" => rho_hi_ts, + "T_probe" => T_probe, + "beta" => BETA, + "Kp" => KP_HEATUP, + "T_c0" => T_C0, + "T_cold0" => T_COLD0, + "T_standby" => T_STANDBY, +)) +println("\nSaved full per-step envelopes to $mat_out") diff --git a/docs/figures/reach_operation_tubes.png b/docs/figures/reach_operation_tubes.png new file mode 100644 index 0000000..6b23011 Binary files /dev/null and b/docs/figures/reach_operation_tubes.png differ diff --git a/reachability/predicates.json b/reachability/predicates.json index b0c8f82..3dff282 100644 --- a/reachability/predicates.json +++ b/reachability/predicates.json @@ -124,6 +124,14 @@ "halfspaces": [ { "row": [[8, -0.4587], [9, 0.9587], [10, -0.5000]], "rhs_expr": "0.01389" } ] + }, + "prompt_critical_margin_heatup": { + "meaning": "Total reactivity stays safely below prompt-critical; keeps the prompt-jump PJ reduction valid. Per-mode specialization: for heatup with feedback linearization, rho_total = Kp*(T_ref(t) - T_c). Requiring rho <= 0.5*beta gives T_c >= T_ref - (beta - 0.5*beta)/Kp = T_ref - 32.5. Since T_ref increases monotonically from T_standby up to T_c0, the worst case is at the ramp end: T_c >= T_c0 - 32.5 = 275.85 C. Approximate as a constant halfspace with the worst-case T_ref.", + "concretization": "T_c >= 275.85 (ensures rho <= 0.5*beta under heatup FL controller)", + "halfspaces": [ + { "state_index": 9, "coeff": -1.0, "rhs_expr": "-275.85" } + ], + "_note": "This is a controller-specific specialization. If the controller or its gain change, recompute. For LQR operation and constant-u scram the margin is trivially satisfied (much larger)." } }, @@ -138,9 +146,10 @@ "fuel_centerline", "cold_leg_subcooled", "heatup_rate_upper", - "heatup_rate_lower" + "heatup_rate_lower", + "prompt_critical_margin_heatup" ], - "_note": "dT_c/dt is linear in (T_f, T_c, T_cold), so ramp-rate halfspace has 3 nonzero coefficients per row. DNBR not modeled — would need an augmented correlation-based predicate." + "_note": "dT_c/dt is linear in (T_f, T_c, T_cold), so ramp-rate halfspace has 3 nonzero coefficients per row. prompt_critical_margin_heatup proves the PJ reduction's validity condition (beta - rho > 0 with margin) as part of the safety obligation rather than as an unverified premise. DNBR not modeled — would need an augmented correlation-based predicate." }, "inv2_holds": { "meaning": "Operation mode safety envelope.",