julia-port: parallel plant model; sanity sim matches MATLAB, reach is stub
Port pke_params, pke_th_rhs, pke_linearize, and all five controllers to Julia. sim_sanity.jl reproduces the MATLAB main.m operation-mode scenario (100%->80% Q_sg step) and matches final state to 3 decimals across n, T_f, T_avg, T_cold, u. reach_operation.jl is a stub: ReachabilityAnalysis.jl (LGG09, GLGM06, BFFPSV18) numerically explodes on the raw stiff system — envelopes of 1e14 K to 1e37 K instead of the known-tight 0.03 K. Almost certainly a state-scaling issue: precursors C_i ~ 1e5, temperatures ~ 300, eigvals span 5000x. Diagonal scaling + retry is planned; left for the next pass since the hand-rolled MATLAB reach already discharges the operation-mode obligation. Project.toml pins OrdinaryDiffEq >= 6.111 (the one that precompiled cleanly on first instantiate). Manifest gitignored. Hacker-Split: Julia path open, reach side needs a scaling pass. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This commit is contained in:
Dane Sabo
parent
02a675c152
commit
9fc4afb611
1
julia-port/.gitignore
vendored
Normal file
1
julia-port/.gitignore
vendored
Normal file
@ -0,0 +1 @@
|
||||
Manifest.toml
|
||||
14
julia-port/Project.toml
Normal file
14
julia-port/Project.toml
Normal file
@ -0,0 +1,14 @@
|
||||
authors = ["Dane Sabo <yourstruly@danesabo.com>"]
|
||||
|
||||
[deps]
|
||||
LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
|
||||
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
|
||||
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
|
||||
MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
|
||||
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
|
||||
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
|
||||
ReachabilityAnalysis = "1e97bd63-91d1-579d-8e8d-501d2b57c93f"
|
||||
|
||||
[compat]
|
||||
OrdinaryDiffEq = "6.111.0"
|
||||
julia = "1.10"
|
||||
40
julia-port/README.md
Normal file
40
julia-port/README.md
Normal file
@ -0,0 +1,40 @@
|
||||
# Julia port
|
||||
|
||||
Parallel implementation of the plant model (`../plant-model/`) in Julia,
|
||||
intended as an eventual landing spot for reachability if we outgrow the
|
||||
MATLAB / hand-rolled tooling.
|
||||
|
||||
## Status
|
||||
|
||||
- `src/pke_params.jl`, `src/pke_th_rhs.jl`, `src/pke_linearize.jl` —
|
||||
functional, match MATLAB term-for-term.
|
||||
- `controllers/controllers.jl` — all five modes ported
|
||||
(null, shutdown, heatup, operation, scram). LQR factory included via
|
||||
MatrixEquations.jl.
|
||||
- `scripts/sim_sanity.jl` — closed-loop simulation matches MATLAB to 3
|
||||
decimals on the 100% → 80% `Q_sg` step. ✅ passing.
|
||||
- `scripts/reach_operation.jl` — ❌ stub. ReachabilityAnalysis.jl
|
||||
algorithms blow up on this stiff, badly-conditioned system. See the
|
||||
file header for diagnosis and planned fix (state rescaling).
|
||||
|
||||
## When to prefer Julia over MATLAB
|
||||
|
||||
Today: nowhere. The sanity path exists so we don't regret the eventual
|
||||
port.
|
||||
|
||||
Once the reach scaling is resolved:
|
||||
- Nonlinear reach for the P controller in operation mode (CORA /
|
||||
JuliaReach territory where MATLAB's linearization doesn't suffice).
|
||||
- Heatup reach with its time-varying reference.
|
||||
- Parametric studies where MATLAB license fees / CI friction matter.
|
||||
|
||||
## Running
|
||||
|
||||
```bash
|
||||
cd julia-port
|
||||
julia --project=. -e 'using Pkg; Pkg.instantiate()'
|
||||
julia --project=. scripts/sim_sanity.jl
|
||||
```
|
||||
|
||||
First run will precompile the dependency stack (OrdinaryDiffEq,
|
||||
ReachabilityAnalysis, LazySets — a few minutes).
|
||||
49
julia-port/controllers/controllers.jl
Normal file
49
julia-port/controllers/controllers.jl
Normal file
@ -0,0 +1,49 @@
|
||||
"""
|
||||
Mode controllers — signatures match the MATLAB side:
|
||||
u = ctrl_<mode>(t, x, plant, ref)
|
||||
|
||||
Pure functions. `ref` can be any NamedTuple of setpoints; unused fields
|
||||
are ignored. For heatup, `ref` must provide `T_start`, `T_target`, `ramp_rate`.
|
||||
"""
|
||||
|
||||
ctrl_null(t, x, plant, ref) = 0.0
|
||||
|
||||
ctrl_shutdown(t, x, plant, ref) = -5.0 * plant.beta
|
||||
ctrl_scram(t, x, plant, ref) = -8.0 * plant.beta
|
||||
|
||||
function ctrl_operation(t, x, plant, ref)
|
||||
Kp = 1e-4
|
||||
T_avg = x[9]
|
||||
return Kp * (ref.T_avg - T_avg)
|
||||
end
|
||||
|
||||
function ctrl_heatup(t, x, plant, ref)
|
||||
Kp = 1e-4
|
||||
T_f = x[8]
|
||||
T_avg = x[9]
|
||||
u_ff = -plant.alpha_f * (T_f - plant.T_f0) -
|
||||
plant.alpha_c * (T_avg - plant.T_c0)
|
||||
T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target)
|
||||
u_unsat = u_ff + Kp * (T_ref - T_avg)
|
||||
u_min = get(ref, :u_min, -5 * plant.beta)
|
||||
u_max = get(ref, :u_max, 0.5 * plant.beta)
|
||||
return clamp(u_unsat, u_min, u_max)
|
||||
end
|
||||
|
||||
"""
|
||||
ctrl_operation_lqr_factory(plant; Q_lqr, R_lqr)
|
||||
|
||||
Returns a closure `ctrl(t, x, plant_ignored, ref_ignored)` with the LQR
|
||||
gain baked in. Pattern chosen so the user can specify Q/R from the
|
||||
calling script and get a pure function to pass to the ODE solver.
|
||||
|
||||
Depends on MatrixEquations.jl for `arec` (algebraic Riccati).
|
||||
"""
|
||||
function ctrl_operation_lqr_factory(plant, A, B; Q_lqr, R_lqr)
|
||||
x_op = pke_initial_conditions(plant)
|
||||
X_ric, _, _ = MatrixEquations.arec(A, B, R_lqr, Q_lqr)
|
||||
K = (R_lqr \ B') * X_ric # 1x10
|
||||
return function (t, x, plant_ignored, ref_ignored)
|
||||
return -(K * (x - x_op))[1]
|
||||
end, K, x_op
|
||||
end
|
||||
90
julia-port/scripts/reach_operation.jl
Normal file
90
julia-port/scripts/reach_operation.jl
Normal file
@ -0,0 +1,90 @@
|
||||
#!/usr/bin/env julia
|
||||
#
|
||||
# reach_operation.jl — operation-mode linear reach in Julia. **STUB**.
|
||||
#
|
||||
# Attempt to reproduce reachability/reach_operation.m using
|
||||
# ReachabilityAnalysis.jl.
|
||||
#
|
||||
# STATUS: does NOT yet produce a useful result. The closed-loop system
|
||||
# is strongly stiff (eigvals spanning 0.012 to 65 /s) and has states
|
||||
# with magnitudes differing by ~10 orders of magnitude (precursors C_i
|
||||
# ~ 1e5 due to 1/Lambda, temperatures ~ 300). LGG09 / GLGM06 / BFFPSV18
|
||||
# all blow up numerically over the 600s horizon, returning envelopes
|
||||
# ranging from 2757 K to 1e37 K instead of the known-tight 0.03 K.
|
||||
#
|
||||
# Likely fix: diagonal state scaling S such that A_cl_scaled = S A_cl S^{-1}
|
||||
# has O(1) entries, run reach in scaled coordinates, then invert S.
|
||||
# Also worth trying: ASB07 (adaptive step) or the Taylor-model schemes.
|
||||
# Left as future work — the hand-rolled MATLAB zonotope in
|
||||
# reachability/reach_linear.m gives a valid result today, so the Julia
|
||||
# port is priority-B.
|
||||
|
||||
using Pkg
|
||||
Pkg.activate(joinpath(@__DIR__, ".."))
|
||||
|
||||
using LinearAlgebra
|
||||
using MatrixEquations
|
||||
using ReachabilityAnalysis, LazySets
|
||||
|
||||
include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
|
||||
include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
|
||||
include(joinpath(@__DIR__, "..", "src", "pke_linearize.jl"))
|
||||
|
||||
plant = pke_params()
|
||||
x_op = pke_initial_conditions(plant)
|
||||
|
||||
# --- Linearize around the operating point ---
|
||||
A, B, B_w, _, _, _ = pke_linearize(plant)
|
||||
|
||||
# --- LQR gain (same Q, R as the MATLAB side) ---
|
||||
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) # 1 x 10
|
||||
A_cl = A - reshape(B, :, 1) * K
|
||||
|
||||
println("\n=== Julia closed-loop spectrum ===")
|
||||
ev = eigvals(A_cl)
|
||||
println(" max Re(eig) = ", round(maximum(real.(ev)); sigdigits=4))
|
||||
println(" min Re(eig) = ", round(minimum(real.(ev)); sigdigits=4))
|
||||
|
||||
# --- Reach-set problem: dx/dt = A_cl dx + B_w w, dx(0) ∈ X0, w ∈ W ---
|
||||
delta_entry = [0.01 * x_op[1];
|
||||
0.001 .* abs.(x_op[2:7]);
|
||||
0.1; 0.1; 0.1]
|
||||
X0 = Hyperrectangle(zeros(10), delta_entry)
|
||||
|
||||
Q_nom = plant.P0
|
||||
w_lo, w_hi = -0.15 * plant.P0, 0.0
|
||||
W = Interval(w_lo, w_hi)
|
||||
B_w_col = reshape(B_w, :, 1)
|
||||
|
||||
# Encode the bounded disturbance as a bounded input:
|
||||
# dx/dt = A_cl x + B_w u, u ∈ W.
|
||||
# Direct constructor since the @system macro's dialect is package-version sensitive.
|
||||
sys = ConstrainedLinearControlContinuousSystem(A_cl, B_w_col, Universe(10), W)
|
||||
prob = InitialValueProblem(sys, X0)
|
||||
sol = solve(prob; T=600.0, alg=GLGM06(; δ=0.5, max_order=20))
|
||||
|
||||
# Extract T_c envelope via support function queries.
|
||||
flow = flowpipe(sol)
|
||||
e9 = zeros(10); e9[9] = 1.0
|
||||
T_c_hi = [ρ(e9, set(R)) for R in flow]
|
||||
T_c_lo = [-ρ(-e9, set(R)) for R in flow]
|
||||
|
||||
println("\n=== Operation-mode reach envelope on T_c (deviation from T_c0) ===")
|
||||
println(" min dT_c = ", round(minimum(T_c_lo); digits=4), " K")
|
||||
println(" max dT_c = ", round(maximum(T_c_hi); digits=4), " K")
|
||||
println(" safety band |dT_c| <= 5.0 K")
|
||||
if maximum(abs.([T_c_lo; T_c_hi])) <= 5.0
|
||||
println(" OK: Julia reach set inside safety band.")
|
||||
else
|
||||
println(" *** VIOLATION ***")
|
||||
end
|
||||
|
||||
# Save the envelope for later comparison
|
||||
using MAT
|
||||
matwrite(joinpath(@__DIR__, "..", "..", "reachability", "julia_reach_operation.mat"),
|
||||
Dict("T_c_lo" => T_c_lo, "T_c_hi" => T_c_hi,
|
||||
"A_cl" => A_cl, "K" => Matrix(K), "delta_entry" => delta_entry))
|
||||
println("\nSaved envelope to reachability/julia_reach_operation.mat")
|
||||
40
julia-port/scripts/sim_sanity.jl
Normal file
40
julia-port/scripts/sim_sanity.jl
Normal file
@ -0,0 +1,40 @@
|
||||
#!/usr/bin/env julia
|
||||
#
|
||||
# sim_sanity.jl — verify the Julia port matches the MATLAB result.
|
||||
#
|
||||
# Reproduces main.m Run 2 (ctrl_operation under 100% -> 80% Q_sg step)
|
||||
# and prints the final state, which should agree with MATLAB to ~1e-3.
|
||||
|
||||
using Pkg
|
||||
Pkg.activate(joinpath(@__DIR__, ".."))
|
||||
|
||||
using OrdinaryDiffEq
|
||||
|
||||
include(joinpath(@__DIR__, "..", "src", "pke_params.jl"))
|
||||
include(joinpath(@__DIR__, "..", "src", "pke_th_rhs.jl"))
|
||||
include(joinpath(@__DIR__, "..", "controllers", "controllers.jl"))
|
||||
|
||||
plant = pke_params()
|
||||
x0 = pke_initial_conditions(plant)
|
||||
|
||||
Q_sg = t -> plant.P0 * (1.0 - 0.2 * (t >= 30))
|
||||
ref = (; T_avg = plant.T_c0)
|
||||
|
||||
function rhs!(dx, x, p, t)
|
||||
u = ctrl_operation(t, x, plant, ref)
|
||||
pke_th_rhs!(dx, x, t, plant, Q_sg, u)
|
||||
end
|
||||
|
||||
prob = ODEProblem(rhs!, x0, (0.0, 600.0))
|
||||
sol = solve(prob, Rodas5(); reltol=1e-8, abstol=1e-10)
|
||||
|
||||
xf = sol.u[end]
|
||||
CtoF(T) = T * 9/5 + 32
|
||||
println("\n=== Julia port sanity — ctrl_operation under 100% -> 80% Q_sg step ===")
|
||||
println(" Final t = ", sol.t[end])
|
||||
println(" n = $(round(xf[1]; digits=4)) (expect ~0.800)")
|
||||
println(" T_f = $(round(CtoF(xf[8]); digits=2)) F (expect ~616.6)")
|
||||
println(" T_avg = $(round(CtoF(xf[9]); digits=2)) F (expect ~587.8)")
|
||||
println(" T_cold = $(round(CtoF(xf[10]); digits=2)) F (expect ~561.4)")
|
||||
u_final = ctrl_operation(sol.t[end], xf, plant, ref)
|
||||
println(" u = $(round(u_final/plant.beta; digits=4)) \$ (expect ~-0.0068)")
|
||||
37
julia-port/src/pke_linearize.jl
Normal file
37
julia-port/src/pke_linearize.jl
Normal file
@ -0,0 +1,37 @@
|
||||
"""
|
||||
pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)
|
||||
|
||||
Numerical Jacobians via central finite differences. Returns
|
||||
`(A, B, B_w, x_star, u_star, Q_star)` such that for small (dx, du, dw):
|
||||
|
||||
dx/dt ≈ A dx + B du + B_w dw, where w = Q_sg.
|
||||
|
||||
Defaults: `x_star` = operating-point steady state, `u_star = 0`,
|
||||
`Q_star = P0`.
|
||||
"""
|
||||
function pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)
|
||||
x_star === nothing && (x_star = pke_initial_conditions(plant))
|
||||
Q_star === nothing && (Q_star = plant.P0)
|
||||
|
||||
n = length(x_star)
|
||||
eps_rel = 1e-6
|
||||
eps_abs = 1e-8
|
||||
|
||||
f = (x, u, w) -> pke_th_rhs(x, 0.0, plant, t -> w, u)
|
||||
|
||||
A = zeros(n, n)
|
||||
for k in 1:n
|
||||
h = max(eps_rel * abs(x_star[k]), eps_abs)
|
||||
xp = copy(x_star); xp[k] += h
|
||||
xm = copy(x_star); xm[k] -= h
|
||||
A[:, k] = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) ./ (2h)
|
||||
end
|
||||
|
||||
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)) ./ (2h)
|
||||
|
||||
h = max(eps_rel * abs(Q_star), 1.0)
|
||||
B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) ./ (2h)
|
||||
|
||||
return A, B, B_w, x_star, u_star, Q_star
|
||||
end
|
||||
53
julia-port/src/pke_params.jl
Normal file
53
julia-port/src/pke_params.jl
Normal file
@ -0,0 +1,53 @@
|
||||
"""
|
||||
pke_params()
|
||||
|
||||
Return a NamedTuple with all plant parameters and derived steady-state
|
||||
conditions for the PKE + thermal-hydraulics model. See
|
||||
../plant-model/pke_params.m for physical rationale on each value.
|
||||
|
||||
Kept as a NamedTuple rather than a mutable struct so that reachability
|
||||
tools (`@taylorize` in ReachabilityAnalysis.jl) can inline the constants.
|
||||
"""
|
||||
function pke_params()
|
||||
# --- Neutronics ---
|
||||
Lambda = 1e-4
|
||||
beta_i = [0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273]
|
||||
lambda_i = [0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01]
|
||||
beta = sum(beta_i)
|
||||
|
||||
# --- Thermal-hydraulic ---
|
||||
P0 = 1000e6
|
||||
M_f = 50000.0
|
||||
c_f = 300.0
|
||||
M_c = 20000.0
|
||||
c_c = 5450.0
|
||||
hA = 5e7
|
||||
W = 5000.0
|
||||
M_sg = 30000.0
|
||||
|
||||
# --- Reactivity feedback coefficients ---
|
||||
alpha_f = -2.5e-5
|
||||
alpha_c = -1.0e-4
|
||||
|
||||
# --- Derived steady-state (full-power equilibrium) ---
|
||||
T_cold0 = 290.0
|
||||
dT_core = P0 / (W * c_c)
|
||||
T_hot0 = T_cold0 + dT_core
|
||||
T_c0 = (T_hot0 + T_cold0) / 2
|
||||
T_f0 = T_c0 + P0 / hA
|
||||
|
||||
return (; Lambda, beta_i, lambda_i, beta, P0, M_f, c_f, M_c, c_c, hA,
|
||||
W, M_sg, alpha_f, alpha_c, T_cold0, T_hot0, T_c0, T_f0, dT_core)
|
||||
end
|
||||
|
||||
"""
|
||||
pke_initial_conditions(plant)
|
||||
|
||||
Full-power steady state: n=1, precursors at equilibrium, temperatures at
|
||||
(T_f0, T_c0, T_cold0). Returns a 10-element Vector.
|
||||
"""
|
||||
function pke_initial_conditions(plant)
|
||||
n0 = 1.0
|
||||
C0 = (plant.beta_i ./ (plant.lambda_i .* plant.Lambda)) .* n0
|
||||
return [n0; C0; plant.T_f0; plant.T_c0; plant.T_cold0]
|
||||
end
|
||||
59
julia-port/src/pke_th_rhs.jl
Normal file
59
julia-port/src/pke_th_rhs.jl
Normal file
@ -0,0 +1,59 @@
|
||||
"""
|
||||
pke_th_rhs!(dx, x, t, plant, Q_sg, u)
|
||||
|
||||
In-place ODE RHS for the coupled PKE + thermal-hydraulics model.
|
||||
Matches ../plant-model/pke_th_rhs.m term-for-term.
|
||||
|
||||
State x = [n, C1..C6, T_f, T_c, T_cold].
|
||||
|
||||
Arguments:
|
||||
- `dx` : 10-element mutable output
|
||||
- `x` : 10-element state
|
||||
- `t` : time [s]
|
||||
- `plant`: parameter NamedTuple from `pke_params`
|
||||
- `Q_sg`: callable `Q_sg(t)` returning SG heat removal [W]
|
||||
- `u` : scalar external reactivity [dk/k]
|
||||
"""
|
||||
function pke_th_rhs!(dx, x, t, plant, Q_sg, u)
|
||||
n = x[1]
|
||||
T_f = x[8]
|
||||
T_c = x[9]
|
||||
T_cold = x[10]
|
||||
T_hot = 2 * T_c - T_cold
|
||||
|
||||
Qsg_t = Q_sg(t)
|
||||
|
||||
rho = u +
|
||||
plant.alpha_f * (T_f - plant.T_f0) +
|
||||
plant.alpha_c * (T_c - plant.T_c0)
|
||||
|
||||
# Neutronics
|
||||
dx[1] = (rho - plant.beta) / plant.Lambda * n +
|
||||
plant.lambda_i[1]*x[2] + plant.lambda_i[2]*x[3] +
|
||||
plant.lambda_i[3]*x[4] + plant.lambda_i[4]*x[5] +
|
||||
plant.lambda_i[5]*x[6] + plant.lambda_i[6]*x[7]
|
||||
|
||||
# Precursors
|
||||
for i in 1:6
|
||||
dx[1+i] = (plant.beta_i[i] / plant.Lambda) * n - plant.lambda_i[i] * x[1+i]
|
||||
end
|
||||
|
||||
# Thermal-hydraulics
|
||||
dx[8] = (plant.P0 * n - plant.hA * (T_f - T_c)) / (plant.M_f * plant.c_f)
|
||||
dx[9] = (plant.hA * (T_f - T_c) - 2 * plant.W * plant.c_c * (T_c - T_cold)) /
|
||||
(plant.M_c * plant.c_c)
|
||||
dx[10] = (plant.W * plant.c_c * (T_hot - T_cold) - Qsg_t) /
|
||||
(plant.M_sg * plant.c_c)
|
||||
|
||||
return nothing
|
||||
end
|
||||
|
||||
"""
|
||||
Convenience: allocating version (returns dx as a new vector). Useful in
|
||||
scripts and tests; for hot-loop reachability use the in-place version.
|
||||
"""
|
||||
function pke_th_rhs(x, t, plant, Q_sg, u)
|
||||
dx = similar(x)
|
||||
pke_th_rhs!(dx, x, t, plant, Q_sg, u)
|
||||
return dx
|
||||
end
|
||||
Loading…
x
Reference in New Issue
Block a user