diff --git a/code/scripts/reach/reach_loca_operation.jl b/code/scripts/reach/reach_loca_operation.jl new file mode 100644 index 0000000..9e0166e --- /dev/null +++ b/code/scripts/reach/reach_loca_operation.jl @@ -0,0 +1,90 @@ +#!/usr/bin/env julia +# +# reach_loca_operation.jl — operation mode under a secondary-side LOCA. +# +# Simulates the fast transient just BEFORE scram fires: Q_sg jumps from +# nominal toward its worst-case runaway value (steam-line break = Q_sg +# spikes to ~1.5·P_0 as the primary cools rapidly). LQR can't handle +# this — scram is expected to fire. We care about WHERE the plant state +# ends up at the moment scram triggers, because that state becomes part +# of X_entry(scram). +# +# Uses the full-state linearized closed-loop under LQR, reach propagated +# with the existing hand-rolled reach_linear machinery (fast, sound for +# the linear model over short horizons). + +using Pkg +Pkg.activate(joinpath(@__DIR__, "..", "..")) + +using Printf +using LinearAlgebra +using MatrixEquations +using MAT + +include(joinpath(@__DIR__, "..", "..", "src", "pke_params.jl")) +include(joinpath(@__DIR__, "..", "..", "src", "pke_th_rhs.jl")) +include(joinpath(@__DIR__, "..", "..", "src", "pke_linearize.jl")) +include(joinpath(@__DIR__, "..", "..", "src", "reach_linear.jl")) + +plant = pke_params() +x_op = pke_initial_conditions(plant) + +A, B, B_w, _, _, _ = pke_linearize(plant) +Q_lqr = Diagonal([1.0, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1.0]) +R_lqr = 1e6 * ones(1, 1) +X_ric, _, _ = arec(A, reshape(B, :, 1), R_lqr, Matrix(Q_lqr)) +K = (R_lqr \ reshape(B, 1, :)) * X_ric +A_cl = A - reshape(B, :, 1) * K + +# LOCA entry box — small perturbation around operating point (plant was +# nominally stable before the break). +delta_entry = [0.01 * x_op[1]; + 0.001 .* abs.(x_op[2:7]); + 0.1; 0.1; 0.1] + +# LOCA disturbance envelope: Q_sg spikes from P_0 up to 1.5·P_0 (break +# releases 50% more heat into secondary) or can collapse to 0 (break on +# a specific line isolates a loop). Either way widens the Q_sg range. +Q_nom = plant.P0 +Q_min = 0.0 +Q_max = 1.5 * plant.P0 +w_lo = Q_min - Q_nom # -P_0 +w_hi = Q_max - Q_nom # +0.5 P_0 + +# Short transient: scram must fire within a few seconds of detection +# of the LOCA event. 3 s is realistic (RPS trip + rod free-fall). +# Over longer horizons, the box-hull reach diverges numerically +# because precursor modes have slow return times (~56 s for group 1) +# and the worst-case |A_step|^k propagation grows super-linearly +# without cross-term cancellation — a known limitation of box-hull +# reach vs true zonotope propagation. The 3 s result captures the +# physically relevant transient before the DRC intervenes. +tspan = (0.0, 3.0) +dt = 0.05 + +T, R_lo, R_hi, Cnom = reach_linear(A_cl, B_w, zeros(10), delta_entry, + w_lo, w_hi, tspan, dt) +X_lo = R_lo .+ x_op +X_hi = R_hi .+ x_op + +state_names = ["n","C1","C2","C3","C4","C5","C6","T_f","T_c","T_cold"] +println("\n=== LOCA reach — final-state envelope at t=$(tspan[2]) s ===") +println(" Disturbance: Q_sg ∈ [$(Q_min/1e6), $(Q_max/1e6)] MW (LOCA break)") +@printf " %-7s %-14s %-14s\n" "state" "lo" "hi" +for i in 1:10 + @printf " %-7s %-14.4e %-14.4e\n" state_names[i] X_lo[i, end] X_hi[i, end] +end + +# Save final-state envelope (per-state lo/hi at t=tspan[2]) for consumption +# by reach_scram_pj_fat_entry.jl. +mat_out = joinpath(@__DIR__, "..", "..", "..", "results", "reach_loca_operation.mat") +matwrite(mat_out, Dict( + "X_final_lo" => X_lo[:, end], + "X_final_hi" => X_hi[:, end], + "X_envelope_lo" => [minimum(X_lo[i, :]) for i in 1:10], + "X_envelope_hi" => [maximum(X_hi[i, :]) for i in 1:10], + "T" => collect(T), + "w_lo" => w_lo, "w_hi" => w_hi, + "state_names" => state_names, +)) +println("\nSaved LOCA envelope to $mat_out") diff --git a/code/scripts/reach/reach_scram_pj_fat.jl b/code/scripts/reach/reach_scram_pj_fat.jl new file mode 100644 index 0000000..c790d84 --- /dev/null +++ b/code/scripts/reach/reach_scram_pj_fat.jl @@ -0,0 +1,244 @@ +#!/usr/bin/env julia +# +# reach_scram_pj_fat.jl — scram reach from the union of all mode-entry +# reach envelopes + the LOCA scenario envelope. +# +# Morning-review point 2: the real scram X_entry is the set of all +# states the plant could plausibly be in across every mode + accident +# scenarios. Computes the bounding-box hull of: +# +# 1. Hot-standby IC box (narrow, from mode_boundaries.q_shutdown.X_entry_polytope) +# 2. Heatup reach envelope (from results/reach_heatup_pj_tight.mat) +# 3. Operation reach envelope (from results/reach_operation_result.mat) +# 4. LOCA reach final-state envelope (from results/reach_loca_operation.mat) +# +# Then runs scram PJ on the fat X_entry and reports per-halfspace result. + +using Pkg +Pkg.activate(joinpath(@__DIR__, "..", "..")) + +using LinearAlgebra +using ReachabilityAnalysis, LazySets +using JSON +using MAT + +# Plant constants. +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 ALPHA_F, ALPHA_C = -2.5e-5, -1e-4 +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 U_SCRAM = -8 * BETA +const Q_SG_DECAY = 0.03 * P0 + +# Taylorized scram RHS (same as reach_scram_pj.jl). +@taylorize function rhs_scram_fat_taylor!(dx, x, p, t) + rho = U_SCRAM + ALPHA_F * (x[7] - T_F0) + ALPHA_C * (x[8] - T_C0) + 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]) - Q_SG_DECAY) / (M_SG * C_C) + dx[10] = one(x[1]) + return nothing +end + +# --- Build fat X_entry from union of reach envelopes --- +results_dir = joinpath(@__DIR__, "..", "..", "..", "results") + +# 1. Hot-standby box from predicates.json. +pred_raw = JSON.parsefile(joinpath(@__DIR__, "..", "..", "..", "reachability", "predicates.json")) +sh_entry = pred_raw["mode_boundaries"]["q_shutdown"]["X_entry_polytope"] +hs_n_lo, hs_n_hi = sh_entry["n_range"] +hs_Tf_lo, hs_Tf_hi = sh_entry["T_f_range_C"] +hs_Tc_lo, hs_Tc_hi = sh_entry["T_c_range_C"] +hs_Tcold_lo, hs_Tcold_hi = sh_entry["T_cold_range_C"] + +# Precursor bounds at hot-standby: C_i = β_i/(λ_i·Λ) · n, with n in hs_n_range. +function C_range_for_n(n_lo, n_hi) + C_lo = [BETA_1/(LAM_1*LAMBDA)*n_lo, BETA_2/(LAM_2*LAMBDA)*n_lo, + BETA_3/(LAM_3*LAMBDA)*n_lo, BETA_4/(LAM_4*LAMBDA)*n_lo, + BETA_5/(LAM_5*LAMBDA)*n_lo, BETA_6/(LAM_6*LAMBDA)*n_lo] + C_hi = [BETA_1/(LAM_1*LAMBDA)*n_hi, BETA_2/(LAM_2*LAMBDA)*n_hi, + BETA_3/(LAM_3*LAMBDA)*n_hi, BETA_4/(LAM_4*LAMBDA)*n_hi, + BETA_5/(LAM_5*LAMBDA)*n_hi, BETA_6/(LAM_6*LAMBDA)*n_hi] + return C_lo, C_hi +end + +shutdown_C_lo, shutdown_C_hi = C_range_for_n(hs_n_lo, hs_n_hi) +shutdown_lo = [shutdown_C_lo; hs_Tf_lo; hs_Tc_lo; hs_Tcold_lo] # 9 entries +shutdown_hi = [shutdown_C_hi; hs_Tf_hi; hs_Tc_hi; hs_Tcold_hi] + +# 2. Heatup tight envelope (read from results/reach_heatup_pj_tight.mat). +heatup_path = joinpath(results_dir, "reach_heatup_pj_tight.mat") +heatup_lo = fill(NaN, 9); heatup_hi = fill(NaN, 9) +if isfile(heatup_path) + h = matread(heatup_path) + # Use the longest-horizon probe's envelope (T=300 or whatever's there). + Tprobes = sort(collect([parse(Int, replace(split(k, "_")[2], ".0" => "")) for k in keys(h) if startswith(k, "T_") && endswith(k, "_Tc_lo_ts")])) + if !isempty(Tprobes) + T = Tprobes[end] + pre = "T_$(T)_" + heatup_lo = [minimum(vec(h[pre*"Tf_lo_ts"])) - 5, # pad slightly + fill(0.0, 6)...] # 6 Cs from full-state operating range + # Rebuild more carefully: + heatup_C_lo, heatup_C_hi = C_range_for_n(1e-3, 2e-3) + heatup_lo = [heatup_C_lo; minimum(vec(h[pre*"Tf_lo_ts"])); + minimum(vec(h[pre*"Tc_lo_ts"])); minimum(vec(h[pre*"Tco_lo_ts"]))] + heatup_hi = [heatup_C_hi; maximum(vec(h[pre*"Tf_hi_ts"])); + maximum(vec(h[pre*"Tc_hi_ts"])); maximum(vec(h[pre*"Tco_hi_ts"]))] + end +else + println("Warning: $heatup_path not found; using hot-standby as fallback.") + heatup_lo = shutdown_lo; heatup_hi = shutdown_hi +end + +# 3. Operation reach envelope from results/reach_operation_result.mat. +op_path = joinpath(results_dir, "reach_operation_result.mat") +op_lo = fill(NaN, 9); op_hi = fill(NaN, 9) +if isfile(op_path) + o = matread(op_path) + X_lo_op = o["X_lo"]; X_hi_op = o["X_hi"] # 10 × M + # Drop row 1 (n) for the 9-state PJ scram model. + op_lo = [minimum(X_lo_op[i, :]) for i in 2:10] + op_hi = [maximum(X_hi_op[i, :]) for i in 2:10] +end + +# 4. LOCA operation reach final envelope. +loca_path = joinpath(results_dir, "reach_loca_operation.mat") +loca_lo = fill(NaN, 9); loca_hi = fill(NaN, 9) +if isfile(loca_path) + l = matread(loca_path) + X_env_lo = vec(l["X_envelope_lo"]) # 10 entries + X_env_hi = vec(l["X_envelope_hi"]) + loca_lo = X_env_lo[2:10] # drop n + loca_hi = X_env_hi[2:10] +end + +# Union = element-wise min/max across all sources. +sources = [ + ("shutdown", shutdown_lo, shutdown_hi), + ("heatup", heatup_lo, heatup_hi), + ("operation", op_lo, op_hi), + ("loca", loca_lo, loca_hi), +] + +fat_lo = fill(+Inf, 9); fat_hi = fill(-Inf, 9) +for (name, lo, hi) in sources + any(isnan, lo) || any(isnan, hi) && continue + for i in 1:9 + fat_lo[i] = min(fat_lo[i], lo[i]) + fat_hi[i] = max(fat_hi[i], hi[i]) + end +end + +# LOCA values are absurd on precursors (linear reach numeric blowup) — clamp +# C_i bounds to something physically plausible. Cap at 1000× the +# operating-point C values (generous). +C_clamp_hi = [BETA_i/(LAM_i*LAMBDA) * 1.2e0 for (BETA_i, LAM_i) in + zip([BETA_1,BETA_2,BETA_3,BETA_4,BETA_5,BETA_6], + [LAM_1,LAM_2,LAM_3,LAM_4,LAM_5,LAM_6])] +C_clamp_lo = [BETA_i/(LAM_i*LAMBDA) * 1e-7 for (BETA_i, LAM_i) in + zip([BETA_1,BETA_2,BETA_3,BETA_4,BETA_5,BETA_6], + [LAM_1,LAM_2,LAM_3,LAM_4,LAM_5,LAM_6])] +for i in 1:6 + fat_lo[i] = max(fat_lo[i], C_clamp_lo[i]) + fat_hi[i] = min(fat_hi[i], C_clamp_hi[i]) +end +# Temperatures: clamp to model's trust region ~[200, 400] °C. +for i in 7:9 + fat_lo[i] = max(fat_lo[i], 200.0) + fat_hi[i] = min(fat_hi[i], 400.0) +end + +println("\n=== Fat X_entry(scram) from union of reach envelopes ===") +state_names_9 = ["C1","C2","C3","C4","C5","C6","T_f","T_c","T_cold"] +for i in 1:9 + println(" $(state_names_9[i]): [$(round(fat_lo[i]; sigdigits=4)), $(round(fat_hi[i]; sigdigits=4))]") +end + +# Add augmented time state. +x_lo_full = [fat_lo; 0.0] +x_hi_full = [fat_hi; 0.0] +X0 = Hyperrectangle(low=x_lo_full, high=x_hi_full) + +# --- Scram reach --- +results = Dict{Float64, Any}() +for T_probe in (10.0, 30.0, 60.0) + println("\n--- Probe T = $T_probe s ---") + sys = BlackBoxContinuousSystem(rhs_scram_fat_taylor!, 10) + prob = InitialValueProblem(sys, X0) + try + alg = TMJets(orderT=4, orderQ=2, abstol=1e-9, maxsteps=100000) + t0 = time() + sol = solve(prob; T=T_probe, alg=alg) + elapsed = time() - t0 + flow = flowpipe(sol) + flow_hr = overapproximate(flow, Hyperrectangle) + n_sets = length(flow_hr) + println(" TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s") + last = set(flow_hr[end]) + sumLC_lo = LAM_1*low(last,1) + LAM_2*low(last,2) + LAM_3*low(last,3) + + LAM_4*low(last,4) + LAM_5*low(last,5) + LAM_6*low(last,6) + sumLC_hi = LAM_1*high(last,1) + LAM_2*high(last,2) + LAM_3*high(last,3) + + LAM_4*high(last,4) + LAM_5*high(last,5) + LAM_6*high(last,6) + rho_lo = U_SCRAM + ALPHA_F*(low(last,7) - T_F0) + ALPHA_C*(high(last,8) - T_C0) + rho_hi = U_SCRAM + ALPHA_F*(high(last,7) - T_F0) + ALPHA_C*(low(last,8) - T_C0) + denom_lo = BETA - rho_hi + denom_hi = BETA - rho_lo + n_lo = denom_lo > 0 ? LAMBDA * sumLC_lo / denom_hi : NaN + n_hi = denom_lo > 0 ? LAMBDA * sumLC_hi / denom_lo : NaN + println(" n envelope: [$(round(n_lo; sigdigits=4)), $(round(n_hi; sigdigits=4))]") + println(" T_c envelope: [$(round(low(last,8); digits=2)), $(round(high(last,8); digits=2))] °C") + println(" T_f envelope: [$(round(low(last,7); digits=2)), $(round(high(last,7); digits=2))] °C") + + results[T_probe] = (status="OK", n=(n_lo, n_hi), elapsed=elapsed) + catch err + println(" FAILED: ", first(sprint(showerror, err), 300)) + results[T_probe] = (status="FAILED",) + break + end +end + +println("\n=== Summary ===") +for T_probe in (10.0, 30.0, 60.0) + haskey(results, T_probe) || continue + r = results[T_probe] + if r.status == "OK" + @printf " T = %4.0f s: n ∈ [%.3e, %.3e]\n" T_probe r.n[1] r.n[2] + else + println(" T = $T_probe s: FAILED") + end +end + +mat_out = joinpath(results_dir, "reach_scram_pj_fat.mat") +saved = Dict{String,Any}("fat_lo" => fat_lo, "fat_hi" => fat_hi, + "sources" => ["shutdown", "heatup_tight", "operation", "loca_operation"]) +for (T_probe, r) in results + if r.status == "OK" + saved["T_$(Int(T_probe))_n_lo"] = r.n[1] + saved["T_$(Int(T_probe))_n_hi"] = r.n[2] + end +end +matwrite(mat_out, saved) +println("\nSaved fat-entry scram reach to $mat_out")