fbbaebff9f
Full toolchain port. Numerical equivalence verified against MATLAB:
- main_mode_sweep.jl: every mode's final state matches MATLAB to 3-4 dp
- reach_operation.jl: per-halfspace margins match MATLAB exactly
- barrier_lyapunov.jl: per-halfspace bounds match (best Qbar from sweep
yields max|dT_c| = 33.228 K either side)
- barrier_compare_OL_CL.jl: OL gamma 1.038e13, CL gamma 1.848e4
matching the MATLAB result; LQR helps by ~20,000x on every halfspace.
Phase summary:
Phase 1: pke_solver.jl, plot_pke_results.jl (Plots.jl), main_mode_sweep.jl
Phase 2: reach_linear.jl, reach_operation.jl, barrier_lyapunov.jl,
barrier_compare_OL_CL.jl, load_predicates.jl
Phase 3 (this commit): delete plant-model/ entirely, delete reach
code from reachability/ keeping predicates.json + docs,
git mv julia-port/ -> code/, update root README + CLAUDE,
write code/CLAUDE.md and code/README.md, update reach
README + WALKTHROUGH file paths, journal preamble note
that pre-port entries reference MATLAB paths.
Why now: prompt-neutron stiffness in nonlinear reach made it clear we
need TMJets, which is Julia. Already had the Julia plant model
working and matching MATLAB. Two languages = two sources of truth =
two places to drift. One language, one truth.
Manifest.toml gitignored. .mat results gitignored.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
50 lines
1.5 KiB
Julia
50 lines
1.5 KiB
Julia
"""
|
|
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
|