2bbb1871cc
Architecture restructure from morning review: 1. code/scripts/ subdivided into sim/, reach/, barrier/, plot/. Easier nav; `barrier/` is the natural place for SOS scale-up scripts. 2. Heatup PJ reach variants consolidated behind TOML configs. reach_heatup_pj.jl now takes `--config path/to/config.toml`; configs/heatup/baseline.toml (wide entry, from predicates.json) and configs/heatup/tight.toml (narrow entry, reproduces all-6-halfspaces discharged result). Old reach_heatup_pj_tight.jl and reach_heatup_pj_tight_full.jl deleted (superseded). 3. Reach output .mat files moved from reachability/ to results/. reachability/ now = specs + docs; results/ = ephemeral outputs (gitignored *.mat). README added. 4. OVERNIGHT_NOTES.md archived to claude_memory/2026-04-20-21-overnight- session-summary.md (date range in the filename makes the history clearer). All include() / Pkg.activate() paths in scripts updated for the new depth. Smoke tests pass (reach_operation.jl generates its .mat in the new results/ location; sim_sanity.jl matches MATLAB). Presentation outline for the 20-min prelim talk landed in presentations/prelim-presentation/outline.md. 14-slide assertion- evidence format targeting OT-informed cybersecurity audience. Each slide: one declarative assertion + one figure. Outline includes which figures already exist and which need to be created, timing checkpoints, cybersecurity angle to emphasize, and Q&A prep. New config configs/heatup/with_steam_dump.toml + its companion scripts/reach/reach_heatup_pj_sd.jl (12-state RHS with Q_sg as an augmented bounded parameter x[10] and time as x[11]). Kicks off point 3 from morning review. Next up: scram X_entry expansion (morning point 2) — LOCA scenario + union of mode reach envelopes. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
179 lines
6.9 KiB
Julia
179 lines
6.9 KiB
Julia
#!/usr/bin/env julia
|
||
#
|
||
# reach_heatup_pj_sd.jl — heatup PJ reach WITH bounded Q_sg steam dump.
|
||
#
|
||
# Adds an augmented state x[11] = Q_sg as a bounded parameter (dx[11] = 0).
|
||
# The reach propagates the entry-box range of Q_sg forward; at each reach
|
||
# set, the Q_sg extent is the disturbance envelope the controller has to
|
||
# reject. Physical story: operator controls secondary steam dump; actual
|
||
# value unknown but bounded in [0, 0.05·P_0].
|
||
#
|
||
# Time is carried as x[12] = t (dt/dt = 1).
|
||
#
|
||
# Diagnostic script — uses the tight-entry + steam-dump config.
|
||
|
||
using Pkg
|
||
Pkg.activate(joinpath(@__DIR__, "..", ".."))
|
||
|
||
using LinearAlgebra
|
||
using ReachabilityAnalysis, LazySets
|
||
using JSON
|
||
using MAT
|
||
using TOML
|
||
|
||
# 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 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
|
||
|
||
# 12-state RHS: [C_1..C_6, T_f, T_c, T_cold, Q_sg, t]
|
||
# dx[10] = 0 (Q_sg is a bounded parameter)
|
||
# dx[11] = 1 (time)
|
||
# Hmm — indexing conflict with original version. Reorganize:
|
||
# x[1..6] = C_1..C_6
|
||
# x[7] = T_f
|
||
# x[8] = T_c
|
||
# x[9] = T_cold
|
||
# x[10] = Q_sg (constant-parameter, dx[10]=0)
|
||
# x[11] = t (dx[11]=1)
|
||
@taylorize function rhs_heatup_sd_taylor!(dx, x, p, t)
|
||
rho = KP_HEATUP * (T_STANDBY + RAMP_RATE_CS * x[11] - 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]) - x[10]) / (M_SG * C_C) # Q_sg now enters!
|
||
|
||
dx[10] = zero(x[1]) # Q_sg constant
|
||
dx[11] = one(x[1]) # time
|
||
|
||
return nothing
|
||
end
|
||
|
||
function main()
|
||
config_path = length(ARGS) > 0 ? ARGS[1] :
|
||
joinpath(@__DIR__, "..", "..", "configs", "heatup", "with_steam_dump.toml")
|
||
if !isfile(config_path)
|
||
alt = joinpath(@__DIR__, "..", "..", config_path)
|
||
isfile(alt) && (config_path = alt)
|
||
end
|
||
config = TOML.parsefile(config_path)
|
||
|
||
println("\n=== Heatup PJ reach with steam dump — config: $(config["name"]) ===")
|
||
println(" $(get(config, "description", ""))")
|
||
|
||
e = config["entry"]
|
||
n_lo, n_hi = e["n_range"]
|
||
T_f_lo, T_f_hi = e["T_f_range_C"]
|
||
T_c_lo, T_c_hi = e["T_c_range_C"]
|
||
T_cold_lo, T_cold_hi = e["T_cold_range_C"]
|
||
sd = e["steam_dump"]
|
||
Q_lo = sd["Q_lo_fraction_P0"] * P0
|
||
Q_hi = sd["Q_hi_fraction_P0"] * P0
|
||
|
||
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, Q_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, Q_hi, 0.0]
|
||
|
||
X0 = Hyperrectangle(low=x_lo, high=x_hi)
|
||
println(" n ∈ [$n_lo, $n_hi], T_c ∈ [$T_c_lo, $T_c_hi] °C")
|
||
println(" Q_sg (steam dump) ∈ [$(Q_lo/1e6), $(Q_hi/1e6)] MW")
|
||
|
||
t_cfg = config["tmjets"]
|
||
probes = config["probes"]["horizons_seconds"]
|
||
|
||
results = Dict{Float64, Any}()
|
||
for T_probe in probes
|
||
println("\n--- Probe T = $T_probe s ---")
|
||
sys = BlackBoxContinuousSystem(rhs_heatup_sd_taylor!, 11)
|
||
prob = InitialValueProblem(sys, X0)
|
||
try
|
||
alg = TMJets(orderT=t_cfg["orderT"], orderQ=t_cfg["orderQ"],
|
||
abstol=t_cfg["abstol"], maxsteps=t_cfg["maxsteps"])
|
||
t0 = time()
|
||
sol = solve(prob; T=Float64(T_probe), alg=alg)
|
||
elapsed = time() - t0
|
||
flow_hr = overapproximate(flowpipe(sol), Hyperrectangle)
|
||
n_sets = length(flow_hr)
|
||
println(" TMJets: $n_sets reach-sets in $(round(elapsed; digits=1)) s")
|
||
|
||
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)
|
||
Tco_lo_env = minimum(low(set(R), 9) for R in flow_hr)
|
||
Tco_hi_env = maximum(high(set(R), 9) 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)
|
||
println(" T_c envelope: [$(round(Tc_lo_env; digits=2)), $(round(Tc_hi_env; digits=2))] °C")
|
||
println(" T_f envelope: [$(round(Tf_lo_env; digits=2)), $(round(Tf_hi_env; digits=2))] °C")
|
||
println(" T_cold env: [$(round(Tco_lo_env; digits=2)), $(round(Tco_hi_env; digits=2))] °C")
|
||
|
||
# Compare low-T_avg trip (280 °C).
|
||
low_trip_ok = Tc_lo_env >= 280.0
|
||
println(" Low-T_avg trip (T_c ≥ 280): $(low_trip_ok ? "✅" : "× loose")")
|
||
|
||
results[T_probe] = (status="OK", Tc=(Tc_lo_env, Tc_hi_env),
|
||
Tf=(Tf_lo_env, Tf_hi_env),
|
||
Tcold=(Tco_lo_env, Tco_hi_env),
|
||
low_trip_ok=low_trip_ok)
|
||
catch err
|
||
println(" FAILED: ", first(sprint(showerror, err), 300))
|
||
results[T_probe] = (status="FAILED",)
|
||
break
|
||
end
|
||
end
|
||
|
||
mat_out = joinpath(@__DIR__, "..", "..", "..", "results", config["output"]["result_file"])
|
||
saved = Dict{String, Any}("config_name" => config["name"],
|
||
"Q_lo" => Q_lo, "Q_hi" => Q_hi)
|
||
for (T_probe, r) in results
|
||
if r.status == "OK"
|
||
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]
|
||
end
|
||
end
|
||
matwrite(mat_out, saved)
|
||
println("\nSaved to $mat_out")
|
||
end
|
||
|
||
main()
|