julia migration: port MATLAB to Julia, delete MATLAB, rename julia-port -> code
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>
This commit is contained in:
Dane Sabo
parent
fa45e96fd1
commit
fbbaebff9f
49
CLAUDE.md
49
CLAUDE.md
@ -41,8 +41,9 @@ hybrid controller. This repo composes three layers of the same story:
|
||||
|---|---|---|---|
|
||||
| Story | `thesis/` (submodule) | The PhD proposal: math, prose, claims | Motivates and defines the methodology |
|
||||
| Discrete | `fret-pipeline/` | FRET requirements → LTL → ltlsynt → AIGER circuit → state machine | Implements Thrust 1 + 2 |
|
||||
| Continuous | `plant-model/` | 10-state PKE + thermal-hydraulics PWR model (MATLAB/Octave) | Plant for Thrust 3 (reachability / barrier certs) |
|
||||
| Verification | `reachability/` | Continuous-mode verification artifacts | Thrust 3 output — currently empty |
|
||||
| Continuous + Verification | `code/` + `reachability/` | 10-state PKE + thermal-hydraulics PWR model, controllers, reach machinery (all Julia) plus `predicates.json` and `WALKTHROUGH.md` | Plant + Thrust 3 (reachability / barrier certs) |
|
||||
| Journal | `journal/` | LaTeX invention log, dated entries | Project history / re-derivability |
|
||||
| App | `app/` | Pluto.jl predicate explorer | FRET-adjacent UI for hybrid-systems groups |
|
||||
| Hardware | `hardware/` | Ovation HIL artifacts | Integration milestone — currently empty |
|
||||
| Deliverables | `presentations/` | Conference papers and talks (submodules) | Dissemination |
|
||||
|
||||
@ -57,24 +58,31 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM).
|
||||
`synthesize.sh` → re-trace via `trace_aiger.py`. New `.svg` / `.dot` /
|
||||
`.png` land in `fret-pipeline/diagrams/`. Copy the figure you want to cite
|
||||
into `docs/figures/` so the thesis can pick it up.
|
||||
2. **Plant change.** Edit `plant-model/pke_params.m` or dynamics in
|
||||
`pke_th_rhs.m` → re-run `main.m` to verify behavior. Steady-state values
|
||||
feed the `T_avg`, `P_crit` thresholds in the FRET requirements — if you
|
||||
change them here, check the requirement predicates match.
|
||||
3. **Verification work.** Continuous-mode verification (reachability, barrier
|
||||
certificates) reads the guards from `fret-pipeline/specs/synthesis_config_v3.json`
|
||||
and the dynamics from `plant-model/pke_th_rhs.m`. Results go in `reachability/`.
|
||||
2. **Plant change.** Edit `code/src/pke_params.jl` or dynamics in
|
||||
`code/src/pke_th_rhs.jl` → re-run `code/scripts/main_mode_sweep.jl`
|
||||
to verify behavior. Steady-state values feed the `T_avg`, `P_crit`
|
||||
thresholds via `reachability/predicates.json` — if you change
|
||||
parameters here, the predicate concretization may need to follow.
|
||||
3. **Verification work.** Continuous-mode verification reads the
|
||||
FRET-side predicates from `reachability/predicates.json` and the
|
||||
dynamics from `code/src/pke_th_rhs.jl`. Reach scripts in
|
||||
`code/scripts/reach_*.jl` and `barrier_*.jl`. Results land back in
|
||||
`reachability/*.mat` (gitignored).
|
||||
4. **Thesis edit.** `cd thesis && <edit> && git commit && git push origin main`,
|
||||
then from the umbrella root: `git add thesis && git commit` to bump the
|
||||
submodule pointer.
|
||||
|
||||
## Where to look for what
|
||||
|
||||
- Reactor physics, steady-state, ODE dynamics → `plant-model/`
|
||||
- Reactor physics, steady-state, ODE dynamics, controllers, reach, barrier → `code/`
|
||||
- Predicate concretizations (the seam between FRET and reach) → `reachability/predicates.json`
|
||||
- Reach analysis writeup → `reachability/WALKTHROUGH.md`
|
||||
- Controller synthesis, requirements, state machines → `fret-pipeline/`
|
||||
- Research narrative, math, claims, citations → `thesis/`
|
||||
- Shared figures (for the thesis or talks) → `docs/figures/`
|
||||
- How the layers compose → `docs/architecture.md`
|
||||
- Project history / invention log → `journal/`
|
||||
- Predicate-explorer GUI → `app/`
|
||||
|
||||
## Submodule rituals
|
||||
|
||||
@ -92,10 +100,13 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM).
|
||||
- FRET variables: `control_<group> = q_<value>` for modes; everything else is
|
||||
an environment input. See `fret-pipeline/README.md` for the full naming rule.
|
||||
- Plant model units: internal SI (W, kg, °C); printed/plotted in °F.
|
||||
- Continuous state vector: `y = [n; C1..C6; T_f; T_c; T_cold]`. Guards in
|
||||
FRET should reference the algebraic outputs (`T_hot`, `T_avg`) or direct
|
||||
states, not derived quantities.
|
||||
- Continuous state vector: `x = [n; C1..C6; T_f; T_c; T_cold]`. Guards
|
||||
in FRET should reference the algebraic outputs (`T_hot`, `T_avg`) or
|
||||
direct states, not derived quantities.
|
||||
- `Q_sg(t)` is the bounded disturbance — never treat `T_cold` as an input.
|
||||
- The continuous toolchain is **Julia**. The earlier MATLAB
|
||||
implementation was ported on 2026-04-20 and deleted; recover via
|
||||
git history if archaeologically needed.
|
||||
|
||||
## What's deliberately missing (so you don't go looking)
|
||||
|
||||
@ -105,5 +116,13 @@ Figure 1 (Cold Shutdown → Heatup → Power Operation, with SCRAM).
|
||||
diagram → hand-translate to Stateflow. Roadmapped.
|
||||
- Hardware bring-up (Ovation DCS) — scheduled for the integration thrust.
|
||||
`hardware/` is a placeholder.
|
||||
- Continuous-mode verification — `reachability/` is empty. Choice of tool
|
||||
(CORA / Flow* / SpaceEx / JuliaReach) still pending.
|
||||
- Continuous-mode verification — partial. Linear reach for operation
|
||||
mode discharges all 6 `inv2_holds` halfspaces (approximate, not
|
||||
sound). Lyapunov-ellipsoid barrier fails structurally
|
||||
(anisotropy ceiling — see `reachability/WALKTHROUGH.md`). Nonlinear
|
||||
reach (TMJets) functional to 10 s; needs reduced-order PKE for
|
||||
hours-long horizons.
|
||||
- Polytopic / SOS barriers — flagged in WALKTHROUGH as the next-step
|
||||
alternative to the failing quadratic Lyapunov barrier.
|
||||
- Reduced-order (prompt-jump) PKE — known remedy for nonlinear reach
|
||||
stiffness, not yet implemented.
|
||||
|
||||
85
README.md
85
README.md
@ -1,16 +1,23 @@
|
||||
# pwr-hybrid-3-demo
|
||||
|
||||
Preliminary example for the HAHACS thesis — a verified hybrid controller for
|
||||
a small modular PWR startup. Composes three layers into one demonstrable
|
||||
pipeline:
|
||||
Preliminary example for the HAHACS thesis — a verified hybrid controller
|
||||
for a small modular PWR startup. Composes three layers into one
|
||||
demonstrable pipeline:
|
||||
|
||||
- **Discrete layer** (`fret-pipeline/`): FRET natural-language requirements
|
||||
→ LTL → synthesized AIGER controller → state-machine diagram.
|
||||
- **Continuous layer** (`plant-model/`): 10-state point kinetic equation +
|
||||
thermal-hydraulics PWR model with bounded steam-generator heat removal as
|
||||
the disturbance input.
|
||||
- **Research context** (`thesis/`): the HAHACS PhD proposal that motivates
|
||||
and formalizes the methodology.
|
||||
- **Discrete layer** (`fret-pipeline/`): FRET natural-language
|
||||
requirements → LTL → synthesized AIGER controller → state-machine
|
||||
diagram.
|
||||
- **Continuous layer** (`code/`): 10-state point kinetic equation +
|
||||
thermal-hydraulics PWR model with bounded steam-generator heat
|
||||
removal as the disturbance input. Controllers, linearization, LQR,
|
||||
reach-tube propagator, Lyapunov barrier — all Julia.
|
||||
- **Verification artifacts** (`reachability/`): predicate
|
||||
concretizations (single source of truth in `predicates.json`) and
|
||||
the standalone reach analysis writeup (`WALKTHROUGH.md`).
|
||||
- **Research context** (`thesis/`): the HAHACS PhD proposal.
|
||||
- **Lab journal** (`journal/`): chronological invention log in LaTeX.
|
||||
- **Predicate explorer app** (`app/`): Pluto.jl notebook bridging
|
||||
FRET predicates and continuous-state halfspaces.
|
||||
|
||||
## Layout
|
||||
|
||||
@ -18,14 +25,15 @@ pipeline:
|
||||
pwr-hybrid-3-demo/
|
||||
CLAUDE.md AI-facing context and architecture map
|
||||
docs/
|
||||
architecture.md How the discrete and continuous layers compose
|
||||
architecture.md How the layers compose
|
||||
figures/ Shared figures for thesis + talks
|
||||
fret-pipeline/ FRET → ltlsynt → AIGER → state machine
|
||||
plant-model/ PWR point kinetics + thermal-hydraulics
|
||||
reachability/ Continuous-mode verification (linear-model tube + Lyapunov barrier attempt; see README)
|
||||
julia-port/ Parallel plant-model port + ReachabilityAnalysis.jl scaffold
|
||||
code/ Plant model, controllers, reach (all Julia)
|
||||
reachability/ predicates.json + WALKTHROUGH.md
|
||||
app/ Pluto.jl predicate explorer
|
||||
journal/ LaTeX lab notebook
|
||||
hardware/ Ovation HIL artifacts (TBD)
|
||||
claude_memory/ Session notes by AI agents (distilled up into CLAUDE.md over time)
|
||||
claude_memory/ Short AI-context notes
|
||||
thesis/ [submodule] PhD proposal
|
||||
presentations/
|
||||
2026DICE/ [submodule] DICE 2026 abstract
|
||||
@ -50,37 +58,48 @@ python3 scripts/trace_aiger.py circuits/PWR_HYBRID_3_DRC.aag diagrams
|
||||
dot -Tpng diagrams/PWR_HYBRID_3_DRC_states.dot -o diagrams/PWR_HYBRID_3_DRC_states.png
|
||||
```
|
||||
|
||||
Run the plant model (MATLAB in `plant-model/` — Octave compatibility not tested since the LQR pieces landed):
|
||||
Run the plant model and reach analysis:
|
||||
|
||||
```matlab
|
||||
main % original single-scenario demo (null vs operation)
|
||||
main_mode_sweep % all five DRC modes back-to-back, writes to ../docs/figures/
|
||||
test_linearize % Jacobian sanity check, saves linearization for reach
|
||||
```bash
|
||||
cd code
|
||||
julia --project=. -e 'using Pkg; Pkg.instantiate()' # first time only
|
||||
|
||||
julia --project=. scripts/main_mode_sweep.jl # all 5 DRC modes
|
||||
julia --project=. scripts/reach_operation.jl # operation-mode linear reach
|
||||
julia --project=. scripts/barrier_lyapunov.jl # Lyapunov barrier
|
||||
julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier
|
||||
julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup (10s cap)
|
||||
```
|
||||
|
||||
Run the reach artifacts (`reachability/`):
|
||||
Open the predicate explorer:
|
||||
|
||||
```matlab
|
||||
reach_operation % linear reach tube for operation-mode LQR
|
||||
barrier_lyapunov % Lyapunov-ellipsoid barrier cert attempt (sweeps weights)
|
||||
```bash
|
||||
cd app
|
||||
julia --project=. -e 'using Pluto; Pluto.run()'
|
||||
# Browser opens; navigate to predicate_explorer.jl
|
||||
```
|
||||
|
||||
**Soundness note:** the current reach tube is the LINEAR model's tube;
|
||||
it is not yet a sound over-approximation of the nonlinear plant. See
|
||||
`reachability/README.md` § Soundness status.
|
||||
**Soundness note:** the current reach tubes are over-approximations
|
||||
of the LINEAR model, not sound over-approximations of the nonlinear
|
||||
plant. See `reachability/README.md` and `reachability/WALKTHROUGH.md`.
|
||||
|
||||
## Prerequisites
|
||||
|
||||
- Python 3.10+
|
||||
- [Spot](https://spot.lre.epita.fr/) for `ltlsynt` (`brew install spot`)
|
||||
- [Graphviz](https://graphviz.org/) for `dot` (`brew install graphviz`)
|
||||
- MATLAB or GNU Octave for the plant model
|
||||
- LaTeX (via `latexmk`) for the thesis submodule
|
||||
- Julia 1.10+ (via `juliaup`).
|
||||
- Python 3.10+ (FRET pipeline only).
|
||||
- [Spot](https://spot.lre.epita.fr/) for `ltlsynt` (`brew install spot`).
|
||||
- [Graphviz](https://graphviz.org/) for `dot` (`brew install graphviz`).
|
||||
- LaTeX (via `latexmk`) for the thesis + journal builds.
|
||||
|
||||
## Further reading
|
||||
|
||||
- `CLAUDE.md` — orientation for AI agents working in this repo
|
||||
- `docs/architecture.md` — how the layers compose
|
||||
- `code/CLAUDE.md` — code architecture, conventions, validity range
|
||||
- `code/README.md` — usage and dependencies
|
||||
- `reachability/README.md` — reach scope, soundness status
|
||||
- `reachability/WALKTHROUGH.md` — standalone analysis writeup
|
||||
- `journal/README.md` — journal format conventions
|
||||
- `journal/journal.tex` — the journal itself, dated entries
|
||||
- `thesis/CLAUDE.md` — the thesis project structure
|
||||
- `fret-pipeline/README.md` — FRET naming conventions and pipeline details
|
||||
- `plant-model/README.md` — scenario setup and model equations
|
||||
|
||||
0
julia-port/.gitignore → code/.gitignore
vendored
0
julia-port/.gitignore → code/.gitignore
vendored
185
code/CLAUDE.md
Normal file
185
code/CLAUDE.md
Normal file
@ -0,0 +1,185 @@
|
||||
# CLAUDE.md
|
||||
|
||||
Guidance for Claude Code (and any AI agent) working in this subdirectory.
|
||||
|
||||
See also: the parent repo's `CLAUDE.md` (living-documentation rule —
|
||||
update this file when new knowledge is created), `claude_memory/` for
|
||||
session notes, and `../journal/` for the long-form invention log.
|
||||
|
||||
## What this is
|
||||
|
||||
The full plant-model + reachability toolchain for the PWR_HYBRID_3
|
||||
preliminary example. All Julia — the earlier MATLAB implementation
|
||||
was ported end-to-end on 2026-04-20 and deleted. A 10-state point
|
||||
kinetic equation (PKE) model of a PWR coupled with lumped-parameter
|
||||
thermal-hydraulics, plus mode-specific continuous controllers, linear
|
||||
reach-set machinery, Lyapunov-barrier attempt, and scaffolding for
|
||||
nonlinear reach via TMJets.
|
||||
|
||||
The model is being developed for hybrid-systems reachability analysis.
|
||||
The end goal is nonlinear reach (sound, not linearized) over per-mode
|
||||
reach-avoid obligations, which will require a reduced-order PKE to
|
||||
eliminate prompt-neutron stiffness.
|
||||
|
||||
## Naming convention (thesis-aligned)
|
||||
|
||||
Continuous state and control-theory notation matches the thesis:
|
||||
|
||||
| Symbol | Meaning |
|
||||
|---|---|
|
||||
| `t` | time [s] |
|
||||
| `x` | continuous state vector, `[n; C1..C6; T_f; T_c; T_cold]` |
|
||||
| `plant` | parameter NamedTuple from `pke_params()` |
|
||||
| `u` | control input: external rod reactivity [dk/k] |
|
||||
| `Q_sg` | bounded disturbance: SG heat removal [W], as a callable `Q_sg(t)` |
|
||||
| `ref` | reference/setpoint NamedTuple (e.g. `ref.T_avg`) |
|
||||
|
||||
`T_avg = T_c` (algebraic alias); `T_hot = 2*T_c - T_cold` (linear
|
||||
coolant profile assumption).
|
||||
|
||||
## Running
|
||||
|
||||
From the project root:
|
||||
|
||||
```bash
|
||||
cd code
|
||||
julia --project=. -e 'using Pkg; Pkg.instantiate()' # first time only
|
||||
|
||||
julia --project=. scripts/main_mode_sweep.jl # runs all 5 DRC modes
|
||||
julia --project=. scripts/reach_operation.jl # operation-mode linear reach
|
||||
julia --project=. scripts/barrier_lyapunov.jl # Lyapunov-ellipsoid barrier
|
||||
julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier bounds
|
||||
julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup reach (10s max, stiffness-limited)
|
||||
julia --project=. scripts/sim_sanity.jl # operation-mode nominal sim
|
||||
julia --project=. scripts/sim_heatup.jl # heatup nominal sim
|
||||
```
|
||||
|
||||
Figures save to `../docs/figures/`. Reach results save to
|
||||
`../reachability/*.mat` (gitignored).
|
||||
|
||||
## Architecture
|
||||
|
||||
```
|
||||
code/
|
||||
Project.toml Julia package manifest
|
||||
README.md usage overview
|
||||
src/
|
||||
pke_params.jl plant parameters and derived steady state
|
||||
pke_th_rhs.jl dynamics f(t, x, plant, Q_sg, u)
|
||||
pke_linearize.jl numerical (A, B, B_w) Jacobians
|
||||
pke_solver.jl closed-loop OrdinaryDiffEq driver
|
||||
plot_pke_results.jl 4-panel results plot (Plots.jl)
|
||||
reach_linear.jl hand-rolled box reach propagator
|
||||
load_predicates.jl reads reachability/predicates.json
|
||||
controllers/
|
||||
controllers.jl ctrl_null, shutdown, heatup, operation,
|
||||
operation_lqr factory, scram
|
||||
ctrl_heatup_unsat.jl saturation-disabled variant for reach
|
||||
scripts/
|
||||
main_mode_sweep.jl runs every mode, saves all figures
|
||||
sim_sanity.jl operation-mode Rodas5 sim (MATLAB-parity check)
|
||||
sim_heatup.jl heatup Rodas5 sim
|
||||
reach_operation.jl linear reach tube + inv2_holds check
|
||||
barrier_lyapunov.jl Lyapunov barrier sweep
|
||||
barrier_compare_OL_CL.jl OL vs CL barrier comparison
|
||||
reach_heatup_nonlinear.jl TMJets nonlinear reach (10s horizon cap)
|
||||
```
|
||||
|
||||
## Controller contract
|
||||
|
||||
```julia
|
||||
u = ctrl_mode(t, x, plant, ref) # all five args required
|
||||
```
|
||||
|
||||
Pure function. Returns scalar `u` (external rod reactivity in `dk/k`).
|
||||
All four signature args are required even if unused — lets the solver
|
||||
swap controllers by function reference without caring which one it is.
|
||||
|
||||
## Solver contract
|
||||
|
||||
```julia
|
||||
t, X, U = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan; x0=nothing, verbose=true)
|
||||
```
|
||||
|
||||
Returns `(t::Vector, X::Matrix, U::Vector)` where `X[k, :]` is the
|
||||
state at time `t[k]` and columns of X match the state vector ordering.
|
||||
The underlying integrator is `Rodas5` from OrdinaryDiffEq (stiff-aware,
|
||||
matches the MATLAB `ode15s` behavior we validated against). Optional
|
||||
`x0` argument for starting from a non-equilibrium IC (heatup, shutdown).
|
||||
|
||||
## Key design decisions
|
||||
|
||||
- **`Q_sg` is the disturbance, never an actuator.** Cold leg
|
||||
temperature is a dynamic state resulting from SG heat removal — the
|
||||
reachability formulation needs `Q_sg ∈ [Q_min, Q_max]` as a bounded
|
||||
disturbance, not a control input.
|
||||
- **`u` is explicit.** `pke_th_rhs!` takes `u` directly instead of
|
||||
reading it from the plant struct.
|
||||
- **Controller gains live in the controller file.** Rough-out phase —
|
||||
once we tune more than one mode we'll lift to a config struct.
|
||||
- **`Rodas5` is required** because the system is stiff:
|
||||
`Lambda` (~10⁻⁴ s) vs thermal time constants (~10–100 s).
|
||||
- **LQR gain is cached** in a `Ref` inside `ctrl_operation_lqr_factory`.
|
||||
Reconstruct the factory (not just call it) to retune.
|
||||
- **Plant parameters live in `pke_params()` as a NamedTuple.** For
|
||||
reach analysis (TMJets) we duplicate them as `const` globals in
|
||||
each reach script — `@taylorize` needs compile-time constants, not
|
||||
struct-field accesses. Keep these in sync.
|
||||
|
||||
## Conventions
|
||||
|
||||
- One function per file where natural; group related functions in
|
||||
one file when they're small (the controllers all live in one file).
|
||||
- Internal model uses SI (W, kg, °C); printed/plotted temperatures
|
||||
in °F.
|
||||
- Reference setpoints in SI (so `ref.T_avg` is in °C).
|
||||
- Controllers are pure functions — no persistent state. If a mode
|
||||
needs integral action, augment `x`, don't hide state in globals.
|
||||
(Exception: the LQR factory caches K in a closure, but that's
|
||||
immutable after construction.)
|
||||
|
||||
## Model validity range
|
||||
|
||||
The feedback coefficients `alpha_f`, `alpha_c` are *linear* slopes
|
||||
taken about the full-power operating point `(T_f0, T_c0)`. They become
|
||||
unphysical when extrapolated far from that reference. True cold
|
||||
shutdown (~50 °C everywhere) violates the linear-coefficient
|
||||
assumption.
|
||||
|
||||
**Trust range:** roughly ±50 °C around the operating point.
|
||||
"Shutdown" in this demo means hot standby (~275 °C per
|
||||
`predicates.json` `T_standby`), not cold shutdown. Extending to
|
||||
cold-shutdown semantics requires piecewise-linear or table-lookup
|
||||
temperature coefficients — separate modeling project.
|
||||
|
||||
## For reachability work
|
||||
|
||||
- `pke_th_rhs!(dx, x, t, plant, Q_sg, u)` is the `ẋ = f(x, u)` with
|
||||
`Q_sg` as the disturbance `w`.
|
||||
- Use `pke_linearize(plant; x_star, u_star, Q_star)` for `(A, B, B_w)`
|
||||
at any operating point.
|
||||
- Safe regions come from `../reachability/predicates.json` — three
|
||||
layers: `operational_deadbands` (DRC transitions),
|
||||
`safety_limits` (hard halfspaces), `mode_invariants` (conjunctions
|
||||
used by per-mode reach).
|
||||
- Do NOT try to verify the whole hybrid system at once. Pick one
|
||||
mode, compute its reach set, discharge the per-mode obligation.
|
||||
- Linear reach in `scripts/reach_operation.jl`; nonlinear (TMJets)
|
||||
in `scripts/reach_heatup_nonlinear.jl`.
|
||||
|
||||
## Robustness caveats (idealized in current artifacts)
|
||||
|
||||
- **α_f, α_c are treated as known exactly.** The feedback-linearization
|
||||
in `ctrl_heatup` assumes the controller's α matches the plant's;
|
||||
real reactors have α drift (burnup, boron, xenon). Robust treatment
|
||||
= α as bounded parametric uncertainty.
|
||||
- **Saturation is a hybrid sub-mode.** `ctrl_heatup`'s `clamp(u, ...)`
|
||||
is formally piecewise affine. Current reach treats it as dormant
|
||||
(true for operation/LQR, near-true for the demo heatup trajectory).
|
||||
A rigorous heatup reach has to model saturation regions explicitly.
|
||||
- **Linear-model reach is not sound for the nonlinear plant.** The
|
||||
linear-reach artifacts produce a sound over-approximation of the
|
||||
LINEAR model's reach set, not of the plant's. To upgrade:
|
||||
nonlinear reach directly (TMJets, once the reduced-order PKE
|
||||
eliminates prompt-neutron stiffness), or Taylor-remainder inflation
|
||||
on linear tubes.
|
||||
59
code/README.md
Normal file
59
code/README.md
Normal file
@ -0,0 +1,59 @@
|
||||
# code
|
||||
|
||||
Plant model, controllers, and reach-analysis toolchain for the
|
||||
HAHACS preliminary example. All Julia.
|
||||
|
||||
## What this is
|
||||
|
||||
A 10-state coupled neutronics + thermal-hydraulics model (point kinetic
|
||||
equations + lumped thermal loop) with continuous-mode controllers for
|
||||
each of the DRC modes (shutdown, heatup, operation, scram), plus a
|
||||
hand-rolled linear reach-tube propagator, a Lyapunov-ellipsoid barrier
|
||||
attempt, and scaffolding for TMJets-based nonlinear reach.
|
||||
|
||||
Ported from MATLAB on 2026-04-20 once the reach experiments made it
|
||||
clear that Julia's stack (`OrdinaryDiffEq`, `MatrixEquations`,
|
||||
`ReachabilityAnalysis`, `LazySets`, `@taylorize`) was the right tool
|
||||
for everything we need going forward. The MATLAB originals are in
|
||||
the git history.
|
||||
|
||||
## Running
|
||||
|
||||
First time:
|
||||
|
||||
```bash
|
||||
cd code
|
||||
julia --project=. -e 'using Pkg; Pkg.instantiate()'
|
||||
```
|
||||
|
||||
Subsequent:
|
||||
|
||||
```bash
|
||||
julia --project=. scripts/main_mode_sweep.jl # all 5 DRC modes, figures
|
||||
julia --project=. scripts/reach_operation.jl # operation-mode linear reach
|
||||
julia --project=. scripts/barrier_lyapunov.jl # Lyapunov barrier attempt
|
||||
julia --project=. scripts/barrier_compare_OL_CL.jl # OL vs CL barrier
|
||||
julia --project=. scripts/reach_heatup_nonlinear.jl # nonlinear heatup (10s cap)
|
||||
```
|
||||
|
||||
Figures save to `../docs/figures/`. Reach results save to
|
||||
`../reachability/*.mat` (gitignored).
|
||||
|
||||
## Structure
|
||||
|
||||
See `CLAUDE.md` for the architectural overview and `../journal/` for
|
||||
the invention-log-style narrative of how this code got written.
|
||||
|
||||
## Dependencies
|
||||
|
||||
From `Project.toml`:
|
||||
|
||||
- `OrdinaryDiffEq` — ODE solver, Rodas5 for stiff systems.
|
||||
- `MatrixEquations` — `arec` for LQR Riccati, `lyapc` for Lyapunov.
|
||||
- `ReachabilityAnalysis` + `LazySets` — reach sets and set operations.
|
||||
- `Plots` — figures (GR backend by default).
|
||||
- `JSON` — read `../reachability/predicates.json`.
|
||||
- `MAT` — save results.
|
||||
|
||||
`Manifest.toml` is gitignored; regenerate locally on first
|
||||
`Pkg.instantiate()`.
|
||||
@ -33,6 +33,15 @@ Each section is a dated session. Sessions are written in two styles:
|
||||
content that should be expanded in a later A-pass.
|
||||
\end{itemize}
|
||||
|
||||
\textbf{File-path archaeology note.} Entries dated on or before
|
||||
2026-04-20 (afternoon) refer to MATLAB files under \texttt{plant-model/}
|
||||
and a MATLAB-Julia split where Julia code lived under \texttt{julia-port/}.
|
||||
The 2026-04-20 evening mega-session ported everything to Julia and
|
||||
deleted MATLAB; the result is the unified \texttt{code/} directory.
|
||||
Path references in earlier entries are historically accurate for their
|
||||
date; recover the corresponding sources via \texttt{git log} or the
|
||||
parent commit at the time.
|
||||
|
||||
Callout boxes signal specific content types:
|
||||
|
||||
\begin{derivation}
|
||||
|
||||
@ -1,40 +0,0 @@
|
||||
# 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).
|
||||
@ -1,196 +0,0 @@
|
||||
# CLAUDE.md
|
||||
|
||||
Guidance for Claude Code (and any AI agent) working in this subdirectory.
|
||||
|
||||
See also: the parent repo's `CLAUDE.md` (living-documentation rule — update
|
||||
this file when new knowledge is created) and `claude_memory/` for session
|
||||
notes that haven't yet earned a place here.
|
||||
|
||||
## What this is
|
||||
|
||||
A 10-state point kinetic equation (PKE) model of a PWR coupled with
|
||||
lumped-parameter thermal-hydraulics, plus mode-specific continuous
|
||||
controllers. Implements the *plant* and the *tactical* layer of the HAHACS
|
||||
methodology — see `../thesis/3-research-approach/approach.tex` sections on
|
||||
transitory, stabilizing, and expulsory modes.
|
||||
|
||||
The model is being developed for hybrid-systems reachability analysis. The
|
||||
end goal is barrier-certificate and reachability methods (in `../reachability/`)
|
||||
to prove each continuous mode satisfies its obligations under bounded
|
||||
`Q_sg(t)` variations.
|
||||
|
||||
## Naming convention (thesis-aligned)
|
||||
|
||||
Continuous state and control-theory notation matches the thesis:
|
||||
|
||||
| Symbol | Meaning |
|
||||
|---|---|
|
||||
| `t` | time [s] |
|
||||
| `x` | continuous state vector, `[n; C1..C6; T_f; T_c; T_cold]` |
|
||||
| `plant` | parameter struct from `pke_params()` (was `p` — renamed to avoid collision with `p_i` predicates) |
|
||||
| `u` | control input: external rod reactivity [dk/k] |
|
||||
| `Q_sg` | bounded disturbance: SG heat removal [W], as a function handle `Q_sg(t)` |
|
||||
| `ref` | reference/setpoint struct (e.g. `ref.T_avg`) |
|
||||
|
||||
`T_avg = T_c` (algebraic alias); `T_hot = 2*T_c - T_cold` (linear coolant
|
||||
profile assumption).
|
||||
|
||||
## Running
|
||||
|
||||
From inside MATLAB:
|
||||
|
||||
```matlab
|
||||
main % runs the default closed-loop scenario
|
||||
```
|
||||
|
||||
From a terminal (figures appear as their own windows):
|
||||
|
||||
```bash
|
||||
matlab -nodesktop -nosplash -r "main" # drops at >> prompt; `exit` to quit
|
||||
matlab -r "main" & # full IDE, backgrounded
|
||||
matlab -batch "main" # headless sanity check (no figures shown)
|
||||
```
|
||||
|
||||
`main.m` configures the scenario (plant, `Q_sg`, ref, tspan), picks a
|
||||
controller, and calls `pke_solver`. Tested on MATLAB R2025b; works on any
|
||||
MATLAB with `ode15s` (R2020b or newer).
|
||||
|
||||
## Architecture
|
||||
|
||||
```
|
||||
plant-model/
|
||||
pke_params.m plant parameters and derived steady state
|
||||
pke_initial_conditions.m analytic equilibrium x0 from plant
|
||||
pke_th_rhs.m dynamics: dx/dt = f(t, x, plant, Q_sg, u)
|
||||
pke_linearize.m numerical (A, B, B_w) at any x_star
|
||||
pke_solver.m closed-loop driver (optional x0 arg)
|
||||
plot_pke_results.m 4-panel results plot
|
||||
load_profile.m canned Q_sg demand shapes
|
||||
main.m entry point — single-mode scenarios
|
||||
main_mode_sweep.m runs all DRC modes back-to-back
|
||||
test_linearize.m sanity-check Jacobians vs nonlinear sim
|
||||
controllers/
|
||||
ctrl_null.m u = 0 (baseline, no feedback)
|
||||
ctrl_shutdown.m hot standby: u = -5*beta (constant)
|
||||
ctrl_heatup.m ramp T_avg: feedback lin. + P on ramp, saturated u
|
||||
ctrl_operation.m stabilizing: P on T_avg
|
||||
ctrl_operation_lqr.m full-state LQR around x_op (persistent-cached K)
|
||||
ctrl_scram.m emergency: u = -8*beta (constant)
|
||||
```
|
||||
|
||||
**Controller contract:**
|
||||
|
||||
```matlab
|
||||
u = ctrl_<mode>(t, x, plant, ref)
|
||||
```
|
||||
|
||||
Pure function. Returns scalar `u` (external rod reactivity in `dk/k`). All
|
||||
four signature args are required even if unused — lets the solver swap
|
||||
controllers by function handle without caring which one it is.
|
||||
|
||||
**Solver contract:**
|
||||
|
||||
```matlab
|
||||
[t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan)
|
||||
```
|
||||
|
||||
Builds the closed-loop RHS internally, calls `ode15s`, then reconstructs the
|
||||
control trajectory `U` by re-evaluating `ctrl_fn` at each returned time step.
|
||||
For event-driven or hysteretic controllers, this reconstruction may miss
|
||||
interior solver evaluations — revisit if / when that matters.
|
||||
|
||||
## Key design decisions
|
||||
|
||||
- **`Q_sg` is the disturbance, never an actuator.** Cold leg temperature is
|
||||
a dynamic state resulting from SG heat removal — the reachability
|
||||
formulation needs `Q_sg ∈ [Q_min, Q_max]` as a bounded disturbance, not a
|
||||
control input.
|
||||
- **`u` is explicit.** `pke_th_rhs` takes `u` directly instead of reading
|
||||
`p.rho_ext(t)` from the params struct. This decouples plant from
|
||||
controller and matches the thesis's `f(x, u)` formulation.
|
||||
- **Controller gains live in the controller file.** Rough-out phase — we'll
|
||||
lift to a config struct once we're tuning more than one mode. The gain
|
||||
in `ctrl_operation` (`Kp = 1e-4 /K`) is chosen to roughly double the
|
||||
effective moderator coefficient, not for precision tracking.
|
||||
- **`ode15s` is required** because the system is stiff: `Lambda` (~10⁻⁴ s)
|
||||
vs thermal time constants (~10–100 s).
|
||||
|
||||
## Conventions
|
||||
|
||||
- One function per file, each file has a single responsibility.
|
||||
- Internal model uses SI (W, kg, °C); all printed/plotted temperatures in °F.
|
||||
- Reference setpoints in SI (so `ref.T_avg` is in °C).
|
||||
- Controllers are pure functions — no persistent state. If a mode ever
|
||||
needs integral action, augment `x` rather than hiding state in globals.
|
||||
|
||||
## Model validity range
|
||||
|
||||
The feedback coefficients `alpha_f`, `alpha_c` are *linear* slopes taken
|
||||
about the full-power operating point `(T_f0, T_c0)`. They become
|
||||
unphysical when extrapolated far from that reference — e.g. at true
|
||||
cold shutdown (~50 C everywhere), the model predicts ~+5 $ intrinsic
|
||||
reactivity that real reactors do not exhibit because temperature
|
||||
coefficients saturate and flip sign in the cold-unborated state.
|
||||
|
||||
**Trust range:** roughly ±50 C around the operating point. "Shutdown"
|
||||
in this demo means hot standby (~290 C, n≈0), not cold shutdown.
|
||||
Extending to cold-shutdown semantics requires piecewise-linear or
|
||||
table-lookup temperature coefficients and is out of scope for the
|
||||
running example.
|
||||
|
||||
## For reachability work
|
||||
|
||||
- `pke_th_rhs(t, x, plant, Q_sg, u)` is the `ẋ = f(x, u)` with `Q_sg` as the
|
||||
disturbance `w`. For a given mode, substitute `u = ctrl_<mode>(t, x, plant, ref)`
|
||||
to get the closed-loop `ẋ = f_cl(x, w)`.
|
||||
- Use `pke_linearize(plant, x_star, u_star, Q_star)` for `(A, B, B_w)` at
|
||||
any operating point; returns central-finite-difference Jacobians. The
|
||||
result at the full-power point is saved to
|
||||
`../reachability/linearization_at_op.mat` by `test_linearize.m`.
|
||||
- Safe regions come from the FRET spec's guards (see
|
||||
`../fret-pipeline/specs/synthesis_config_v3.json`) — these define
|
||||
`X_entry`, `X_exit`, `X_safe` per the thesis approach. Numerical
|
||||
predicate thresholds still need to be pinned (see
|
||||
`../claude_memory/` notes for the discussion).
|
||||
- Do NOT try to verify the whole hybrid system at once. Pick one mode,
|
||||
compute its reach set, discharge the per-mode obligation. The
|
||||
compositionality argument in the thesis is what makes this tractable.
|
||||
- First reach set is in `../reachability/reach_operation.m` (operation
|
||||
mode under LQR); a Lyapunov-ellipsoid barrier-cert attempt is in
|
||||
`../reachability/barrier_lyapunov.m`.
|
||||
|
||||
## For control design
|
||||
|
||||
- The LQR gain in `ctrl_operation_lqr.m` is cached in a persistent
|
||||
variable. If you change `Q_lqr` or `R_lqr` or mutate `plant`, restart
|
||||
MATLAB to clear the cache.
|
||||
- `ctrl_heatup.m` uses feedback linearization + saturated P. No
|
||||
integrator on purpose — adding PI would double the plant's intrinsic
|
||||
thermal-mass integration and make anti-windup a hybrid mode, both
|
||||
bad for reach.
|
||||
- Starting heatup from a truly zero-power IC produces a large lag +
|
||||
power-spike overshoot in the simulation. Physical resolution: the
|
||||
reactor should be taken critical in shutdown mode (at ~0.1% power)
|
||||
before DRC transitions to heatup. `main_mode_sweep.m` uses this IC.
|
||||
|
||||
## Robustness caveats (idealized in current artifacts)
|
||||
|
||||
- **α_f, α_c are treated as known exactly.** In reality α_f drifts
|
||||
~20% over burnup; α_c spans ~10x across soluble-boron dilution over
|
||||
a cycle; xenon adds 2-3 $ reactivity on its own timescale. The
|
||||
feedback-linearization in `ctrl_heatup.m` assumes the controller's
|
||||
α matches the plant's; if not, the clean `rho_total = Kp*e`
|
||||
property degrades to `Kp*e + delta*alpha*dT`, and the P term must
|
||||
absorb the residual. Stabilization still holds but reach analysis
|
||||
should eventually treat α as a bounded parametric uncertainty.
|
||||
- **Saturation is a hybrid sub-mode.** The `sat(u, u_min, u_max)` in
|
||||
`ctrl_heatup.m` is formally piecewise affine. Current reach
|
||||
treats it as dormant (true for operation/LQR, near-true for the
|
||||
demo heatup trajectory). A rigorous heatup reach has to model
|
||||
the saturation regions explicitly.
|
||||
- **Linear-model reach is not sound for the nonlinear plant.** The
|
||||
reach artifacts in `../reachability/` use the linearization; the
|
||||
result is a sound over-approximation of the LINEAR model's reach,
|
||||
not of the plant's. To upgrade: nonlinear reach directly, or
|
||||
linear reach + Taylor-remainder inflation. See
|
||||
`../reachability/README.md` § Soundness status.
|
||||
@ -1,92 +0,0 @@
|
||||
# plant-model
|
||||
|
||||
PWR plant model (point kinetics + lumped thermal-hydraulics) and
|
||||
mode-specific continuous controllers for the HAHACS preliminary example.
|
||||
|
||||
## Overview
|
||||
|
||||
A 10-state coupled neutronics + thermal-hydraulics model in MATLAB:
|
||||
|
||||
- 6 delayed neutron precursor groups (U-235 thermal fission, Keepin)
|
||||
- Lumped fuel, core coolant, and SG/cold-leg thermal nodes
|
||||
- Steam generator heat removal `Q_sg(t)` as the bounded disturbance input
|
||||
- Doppler and moderator temperature reactivity feedback
|
||||
- External rod reactivity `u` as the controllable input
|
||||
|
||||
State vector: `x = [n; C1..C6; T_f; T_c; T_cold]` (10 states). See
|
||||
`CLAUDE.md` for the naming convention.
|
||||
|
||||
## Quick Start
|
||||
|
||||
Open MATLAB in this directory and run:
|
||||
|
||||
```matlab
|
||||
main
|
||||
```
|
||||
|
||||
The default scenario runs two simulations of a 100% → 80% SG demand step:
|
||||
once with `ctrl_null` (plant feedback only) and once with `ctrl_operation`
|
||||
(proportional rod reactivity on T_avg error), and plots the comparison.
|
||||
|
||||
## Files
|
||||
|
||||
| File | Role |
|
||||
|------|------|
|
||||
| `main.m` | Entry point — scenario config and run |
|
||||
| `pke_params.m` | Plant parameters and steady-state derivation |
|
||||
| `pke_th_rhs.m` | Dynamics `ẋ = f(t, x, plant, Q_sg, u)` |
|
||||
| `pke_initial_conditions.m` | Analytic steady-state `x0` |
|
||||
| `pke_solver.m` | Closed-loop driver — takes a controller function handle |
|
||||
| `plot_pke_results.m` | 4-panel results plot |
|
||||
| `load_profile.m` | SG heat demand shapes |
|
||||
| `controllers/ctrl_null.m` | `u = 0` baseline |
|
||||
| `controllers/ctrl_operation.m` | Stabilizing mode: P on `T_avg` |
|
||||
|
||||
## Controllers
|
||||
|
||||
Controllers share a single signature:
|
||||
|
||||
```matlab
|
||||
u = ctrl_<mode>(t, x, plant, ref)
|
||||
```
|
||||
|
||||
Returns scalar `u` (external rod reactivity in `dk/k`). The solver swaps
|
||||
controllers via function handle:
|
||||
|
||||
```matlab
|
||||
[t, X, U] = pke_solver(plant, Q_sg, @ctrl_operation, ref, tspan);
|
||||
```
|
||||
|
||||
Additional modes (`ctrl_heatup`, `ctrl_scram`, `ctrl_shutdown`) will land in
|
||||
`controllers/` following the same signature.
|
||||
|
||||
## Running Different Scenarios
|
||||
|
||||
Swap `Q_sg` in `main.m`:
|
||||
|
||||
```matlab
|
||||
% Step down to 90% at t = 10s
|
||||
Q_sg = @(t) plant.P0 * (1.0 - 0.1 * (t >= 10));
|
||||
|
||||
% Interpolated time series
|
||||
t_data = [0, 100, 200, 300];
|
||||
q_data = [1.0, 0.85, 0.9, 1.0] * plant.P0;
|
||||
Q_sg = @(t) interp1(t_data, q_data, t, 'linear', 'extrap');
|
||||
```
|
||||
|
||||
Swap the controller:
|
||||
|
||||
```matlab
|
||||
[t, X, U] = pke_solver(plant, Q_sg, @ctrl_null, [], tspan);
|
||||
```
|
||||
|
||||
Change the reference (for modes that use one):
|
||||
|
||||
```matlab
|
||||
ref.T_avg = plant.T_c0 + 5; % track 5 C above nominal
|
||||
```
|
||||
|
||||
## Requirements
|
||||
|
||||
MATLAB (R2020b or newer, tested on R2025b). Uses `ode15s` from base MATLAB
|
||||
— no toolboxes required.
|
||||
@ -1,88 +0,0 @@
|
||||
# PKE Playground — LLM Context
|
||||
|
||||
This document provides context for an LLM working on this codebase. It explains the model, design decisions, and where the project is headed.
|
||||
|
||||
## What This Is
|
||||
|
||||
A point kinetic equation (PKE) model coupled with lumped-parameter thermal-hydraulics, implemented in MATLAB/Octave. The model represents a simplified PWR (pressurized water reactor) core and primary loop.
|
||||
|
||||
## Why It Exists
|
||||
|
||||
This is a workable plant model being developed for **hybrid systems reachability analysis**. The end goal is to use barrier certificate methods to prove that under bounded variations in steam generator (SG) heat removal, the reactor stays within a defined safe operating region (bounded temperatures, power levels, etc.).
|
||||
|
||||
The model needs to be:
|
||||
- Physically motivated (not just a toy system)
|
||||
- Clean and state-space oriented (explicit state vector, clear inputs/outputs)
|
||||
- Structured so that `Q_sg(t)` is a bounded disturbance input for reachability tools
|
||||
|
||||
## The Model
|
||||
|
||||
### State Vector (10 states)
|
||||
|
||||
```
|
||||
y = [n; C1; C2; C3; C4; C5; C6; T_f; T_c; T_cold]
|
||||
```
|
||||
|
||||
- `n` — normalized neutron population (1.0 = nominal)
|
||||
- `C1..C6` — six delayed neutron precursor group concentrations (Keepin U-235 thermal fission data)
|
||||
- `T_f` — lumped fuel temperature [C]
|
||||
- `T_c` — average core coolant temperature [C]
|
||||
- `T_cold` — cold leg / SG outlet coolant temperature [C]
|
||||
|
||||
### Algebraic Outputs
|
||||
|
||||
- `T_hot = 2*T_c - T_cold` (linear coolant temperature profile through the core)
|
||||
- `T_avg = T_c`
|
||||
|
||||
### Equations
|
||||
|
||||
**Neutronics:**
|
||||
```
|
||||
dn/dt = (rho - beta) / Lambda * n + sum(lambda_i * C_i)
|
||||
dC_i/dt = beta_i / Lambda * n - lambda_i * C_i
|
||||
```
|
||||
|
||||
**Thermal-hydraulics (three lumped nodes):**
|
||||
```
|
||||
M_f*c_f * dT_f/dt = P0*n - hA*(T_f - T_c) [fuel]
|
||||
M_c*c_c * dT_c/dt = hA*(T_f - T_c) - 2*W*c_c*(T_c - T_cold) [core coolant]
|
||||
M_sg*c_c * dT_cold/dt = W*c_c*(T_hot - T_cold) - Q_sg(t) [SG/cold leg]
|
||||
```
|
||||
|
||||
**Reactivity feedback:**
|
||||
```
|
||||
rho = rho_ext(t) + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0)
|
||||
```
|
||||
|
||||
### Why Q_sg Is the Input (Not T_cold)
|
||||
|
||||
In a real PWR, the steam generator removes heat from the primary loop based on secondary-side demand. The cold leg temperature is a *result* of that heat removal, not a boundary condition you get to set. Treating `Q_sg(t)` as the external input and letting `T_cold` be a dynamic state is physically correct and maps naturally to the reachability problem: "given that SG demand is somewhere in [Q_min, Q_max], what states can the reactor reach?"
|
||||
|
||||
### Why T_hot = 2*T_c - T_cold
|
||||
|
||||
This comes from the linear coolant temperature profile assumption through the core. If coolant enters at `T_cold` and exits at `T_hot`, and the average is `T_c`, then `T_c = (T_hot + T_cold) / 2`, so `T_hot = 2*T_c - T_cold`. This avoids needing a separate hot leg state while still tracking both temperatures.
|
||||
|
||||
### Why ode15s
|
||||
|
||||
The PKE system is stiff. The prompt neutron generation time `Lambda` (~10^-4 s) is orders of magnitude smaller than the thermal time constants (~10-100 s) and precursor decay constants (~0.01-3 s^-1). An explicit solver would need extremely small timesteps; `ode15s` handles this naturally.
|
||||
|
||||
## File Structure
|
||||
|
||||
Each file has one purpose:
|
||||
|
||||
| File | Purpose |
|
||||
|------|---------|
|
||||
| `pke_solver.m` | Main driver script — scenario config, solve, print |
|
||||
| `pke_params.m` | All plant parameters and steady-state derivation |
|
||||
| `pke_th_rhs.m` | ODE right-hand side `f(t, y)` — this is the dynamics |
|
||||
| `pke_initial_conditions.m` | Analytic steady-state initial condition |
|
||||
| `load_profile.m` | SG heat demand time series (swappable) |
|
||||
| `plot_pke_results.m` | 4-panel results plotting |
|
||||
|
||||
## For Future Reachability Work
|
||||
|
||||
- `pke_th_rhs.m` is your `dx/dt = f(x, w)` where `w = Q_sg`
|
||||
- Linearize around the steady state from `pke_params.m` to get `A`, `B` matrices
|
||||
- The safe region will be a polytope or sublevel set on `[n, T_f, T_hot, T_cold]` (or a subset)
|
||||
- `Q_sg` is the bounded disturbance: `Q_sg in [Q_min, Q_max]` around `P0`
|
||||
- Barrier certificate `B(x)` would certify that trajectories starting in the safe set stay there for all admissible `Q_sg(t)`
|
||||
@ -1,82 +0,0 @@
|
||||
function u = ctrl_heatup(t, x, plant, ref)
|
||||
% CTRL_HEATUP Ramp T_avg toward a target at a bounded rate.
|
||||
%
|
||||
% Structure:
|
||||
% u_ff = -alpha_f*(T_f - T_f0) - alpha_c*(T_c - T_c0) % cancel feedback
|
||||
% T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target) % ramped reference
|
||||
% u_unsat = u_ff + Kp * (T_ref - T_avg) % P on error
|
||||
% u = sat(u_unsat, ref.u_min, ref.u_max) % bounded rod worth
|
||||
%
|
||||
% Why saturation:
|
||||
% Without it the P gain can push u toward prompt-supercritical as the
|
||||
% cold-hot feedback bias unwinds late in the ramp. Capping u at
|
||||
% +0.5*beta guarantees rho_total < beta (below prompt), which in
|
||||
% turn bounds the neutron-kinetics excursion rate for reachability.
|
||||
%
|
||||
% IMPORTANT for reach analysis:
|
||||
% The sat() is a 3-mode piecewise-affine system (sat-low / linear /
|
||||
% sat-high). Under linear reach assumptions it must either be
|
||||
% (a) proven dormant (u_unsat stays in [u_min, u_max] across the
|
||||
% reach set — trivial to check, expensive to over-approximate
|
||||
% tightly), or
|
||||
% (b) handled as an explicit hybrid automaton nested under the DRC
|
||||
% mode, with transitions when u_unsat crosses the saturation
|
||||
% bounds.
|
||||
% The current reach_operation.m assumes (a) implicitly. Heatup
|
||||
% reach will need option (b) because u_unsat is near the +0.5*beta
|
||||
% bound during early ramp.
|
||||
%
|
||||
% Why no integrator:
|
||||
% Ramp tracking has a structural lag proportional to ramp_rate / Kp_eff.
|
||||
% Acceptable because the DRC exits heatup on a predicate window
|
||||
% (t_avg_in_range & p_above_crit), not on zero steady-state error.
|
||||
% Adding PI would double-count the intrinsic plant integrator
|
||||
% (thermal mass) and make anti-windup a hybrid transition.
|
||||
%
|
||||
% IMPORTANT caveat on the cancellation u_ff:
|
||||
% The feedback linearization works only if the controller's values
|
||||
% of alpha_f, alpha_c match the plant's. In reality alpha drifts:
|
||||
% alpha_f ~20% over burnup, alpha_c by ~10x across boron dilution,
|
||||
% plus xenon. With alpha_true = alpha_nom*(1+delta), the
|
||||
% cancellation leaves residual reactivity delta*alpha*dT that the
|
||||
% P term cleans up — stabilization still holds, but the clean
|
||||
% "rho_total = Kp*e" property is gone. A robust deployment would
|
||||
% treat alpha as an interval and propagate parametric uncertainty
|
||||
% through reach (zonotope with parameter generators), or add
|
||||
% adaptive alpha estimation. Idealized here.
|
||||
%
|
||||
% Inputs:
|
||||
% t - time [s]
|
||||
% x - state vector (10 x 1)
|
||||
% plant - parameter struct (alpha_f, alpha_c, T_f0, T_c0 used)
|
||||
% ref - struct with fields:
|
||||
% .T_start starting T_avg [C]
|
||||
% .T_target final T_avg [C]
|
||||
% .ramp_rate desired dT_avg/dt [C/s]
|
||||
% .u_min (optional) lower saturation [dk/k]; default -5*beta
|
||||
% .u_max (optional) upper saturation [dk/k]; default +0.5*beta
|
||||
|
||||
Kp = 1e-4; % [dk/k per K]
|
||||
|
||||
T_f = x(8);
|
||||
T_c = x(9);
|
||||
T_avg = T_c;
|
||||
|
||||
% Feedforward: cancel intrinsic temperature feedback so rho_total = Kp*e
|
||||
% (before saturation).
|
||||
u_ff = -plant.alpha_f * (T_f - plant.T_f0) ...
|
||||
-plant.alpha_c * (T_avg - plant.T_c0);
|
||||
|
||||
% Ramped reference, clamped at target.
|
||||
T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target);
|
||||
|
||||
e = T_ref - T_avg;
|
||||
u_unsat = u_ff + Kp * e;
|
||||
|
||||
% Saturation bounds (defaults keep rod worth subcritical-prompt).
|
||||
if isfield(ref, 'u_min'), u_min = ref.u_min; else, u_min = -5 * plant.beta; end
|
||||
if isfield(ref, 'u_max'), u_max = ref.u_max; else, u_max = 0.5 * plant.beta; end
|
||||
|
||||
u = min(max(u_unsat, u_min), u_max);
|
||||
|
||||
end
|
||||
@ -1,13 +0,0 @@
|
||||
function u = ctrl_null(t, x, plant, ref)
|
||||
% CTRL_NULL No-op controller. Returns u = 0.
|
||||
%
|
||||
% Use for unforced-plant simulations (feedback-only response to Q_sg) or
|
||||
% as a baseline to compare against closed-loop modes.
|
||||
%
|
||||
% Signature matches all other controllers:
|
||||
% u = ctrl_<mode>(t, x, plant, ref)
|
||||
% so it can be dropped into pke_solver as a function handle.
|
||||
|
||||
u = 0;
|
||||
|
||||
end
|
||||
@ -1,31 +0,0 @@
|
||||
function u = ctrl_operation(t, x, plant, ref)
|
||||
% CTRL_OPERATION Stabilizing mode controller: hold T_avg at ref.T_avg.
|
||||
%
|
||||
% Proportional feedback on external reactivity (rod worth) driven by
|
||||
% temperature error:
|
||||
% u = Kp * (ref.T_avg - T_avg)
|
||||
%
|
||||
% Sign check: when T_avg > ref.T_avg, error < 0, u < 0 → rods insert,
|
||||
% power drops, T_avg falls. Correct polarity for a PWR.
|
||||
%
|
||||
% Gain Kp is a rough-out value. Scale reasoning:
|
||||
% - alpha_c ~ -1e-4 /K means the moderator already supplies 1e-4/K of
|
||||
% restoring reactivity per Kelvin of temperature error.
|
||||
% - Kp = 1e-4 /K roughly doubles the effective moderator coefficient —
|
||||
% it damps transients without overwhelming the intrinsic feedback.
|
||||
% - Real PWR rod controllers use much smaller gains with a deadband.
|
||||
% We can retune or switch to PI once we care about steady-state offset.
|
||||
%
|
||||
% Inputs:
|
||||
% t - time [s] (unused in this controller, kept for signature)
|
||||
% x - state vector (10 x 1)
|
||||
% plant - parameter struct (unused here; reserved for future gain-scheduling)
|
||||
% ref - struct with field .T_avg = target average coolant temperature [C]
|
||||
|
||||
Kp = 1e-4; % [dk/k per K]
|
||||
|
||||
T_avg = x(9);
|
||||
e = ref.T_avg - T_avg;
|
||||
u = Kp * e;
|
||||
|
||||
end
|
||||
@ -1,53 +0,0 @@
|
||||
function u = ctrl_operation_lqr(t, x, plant, ref)
|
||||
% CTRL_OPERATION_LQR Full-state feedback for operation mode.
|
||||
%
|
||||
% u = -K * (x - x_op), K from LQR on (A, B) at x_op.
|
||||
%
|
||||
% Advantages over ctrl_operation (plain P on T_avg):
|
||||
% - Damps the fast neutronics and the slow thermal response jointly.
|
||||
% - Adds NO state (K is 1x10, feedback is memoryless).
|
||||
% - Closed-loop A_cl = A - B*K is trivially propagable for reachability.
|
||||
%
|
||||
% Q/R are persistently cached on the first call. If the caller mutates
|
||||
% plant after the first invocation, gains become stale — clear `persistent`
|
||||
% by restarting MATLAB. That's a reach-analysis bookkeeping issue, not a
|
||||
% runtime one.
|
||||
%
|
||||
% Weighting rationale:
|
||||
% Q diagonal:
|
||||
% n -> 1 (power tracking matters)
|
||||
% C1..C6 -> 1e-3 (precursor deviations are informational)
|
||||
% T_f -> 1e-2 (care about it but not as much as coolant)
|
||||
% T_c -> 1e2 (primary objective)
|
||||
% T_cold -> 1 (secondary but consequential)
|
||||
% R = 1e6 (discourage large rod movements — one Kp of |u| ~ 1e-3)
|
||||
%
|
||||
% Inputs:
|
||||
% t, x, plant, ref - standard signature. ref is unused; the setpoint
|
||||
% is baked in as x_op. If you want to track a
|
||||
% different operating point, regenerate K there.
|
||||
|
||||
persistent K x_op_cached
|
||||
|
||||
if isempty(K)
|
||||
x_op_cached = pke_initial_conditions(plant);
|
||||
[A, B, ~] = pke_linearize(plant, x_op_cached, 0, plant.P0);
|
||||
|
||||
Q = diag([1, ...
|
||||
1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, ...
|
||||
1e-2, 1e2, 1]);
|
||||
R = 1e6;
|
||||
|
||||
% MATLAB's lqr requires Control System Toolbox. Fall back to
|
||||
% icare (in base MATLAB from R2019a) if lqr is unavailable.
|
||||
try
|
||||
K = lqr(A, B, Q, R);
|
||||
catch
|
||||
[~, ~, K_ric] = icare(A, B, Q, R);
|
||||
K = K_ric;
|
||||
end
|
||||
end
|
||||
|
||||
u = -K * (x - x_op_cached);
|
||||
|
||||
end
|
||||
@ -1,20 +0,0 @@
|
||||
function u = ctrl_scram(t, x, plant, ref)
|
||||
% CTRL_SCRAM Emergency shutdown: max negative rod worth inserted.
|
||||
%
|
||||
% Applies the full scram rod worth instantaneously. In a real plant the
|
||||
% rods free-fall in ~2-3 seconds; this idealized version steps to the
|
||||
% final worth at t = t_scram. Good enough for first-pass reachability,
|
||||
% conservative for safety arguments (faster insertion => less excursion).
|
||||
%
|
||||
% u = -8 * beta (~8 $, in the typical 8-10 $ regulatory band)
|
||||
%
|
||||
% Same feedback-linearization comment as ctrl_shutdown: this dominates
|
||||
% any plausible feedback the reactor can produce en route from any
|
||||
% power/temperature state.
|
||||
%
|
||||
% Inputs:
|
||||
% t, x, plant, ref - standard signature; all unused here.
|
||||
|
||||
u = -8 * plant.beta;
|
||||
|
||||
end
|
||||
@ -1,21 +0,0 @@
|
||||
function u = ctrl_shutdown(t, x, plant, ref)
|
||||
% CTRL_SHUTDOWN Hot-standby / cold-shutdown passive controller.
|
||||
%
|
||||
% Holds the reactor deeply subcritical with rods fully in. No feedback —
|
||||
% just a constant negative external reactivity that dominates any
|
||||
% temperature-feedback contribution the reactor could plausibly produce
|
||||
% in its accessible state set.
|
||||
%
|
||||
% u = -5 * beta (~5 $ subcritical)
|
||||
%
|
||||
% Sign check: rho_total = u + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0).
|
||||
% At cold IC (T < T0) the feedback contribution is positive (roughly
|
||||
% +2.8 $ at T = T_cold0 everywhere), so -5 $ leaves the reactor at about
|
||||
% -2.2 $ — comfortably subcritical.
|
||||
%
|
||||
% Inputs:
|
||||
% t, x, plant, ref - standard signature; all unused here.
|
||||
|
||||
u = -5 * plant.beta;
|
||||
|
||||
end
|
||||
@ -1,26 +0,0 @@
|
||||
function f = load_profile(t)
|
||||
% LOAD_PROFILE Returns fractional SG heat demand at time t.
|
||||
%
|
||||
% Output is in [0, 1] where 1.0 = 100% of nominal power P0.
|
||||
% Multiply by P0 to get Q_sg in watts: Q_sg = @(t) P0 * load_profile(t)
|
||||
%
|
||||
% Default profile (load-following):
|
||||
% 0 - 30s: hold at 100%
|
||||
% 30 - 130s: ramp 100% -> 80%
|
||||
% 130 - 350s: hold at 80%
|
||||
% 350 - 450s: ramp 80% -> 100%
|
||||
% 450+: hold at 100%
|
||||
|
||||
if t < 30
|
||||
f = 1.0;
|
||||
elseif t < 130
|
||||
f = 1.0 - 0.2 * (t - 30) / 100;
|
||||
elseif t < 350
|
||||
f = 0.8;
|
||||
elseif t < 450
|
||||
f = 0.8 + 0.2 * (t - 350) / 100;
|
||||
else
|
||||
f = 1.0;
|
||||
end
|
||||
|
||||
end
|
||||
@ -1,50 +0,0 @@
|
||||
%% main.m — closed-loop entry point
|
||||
%
|
||||
% Scenario: a 100% → 80% SG demand step. Compare the plant's intrinsic
|
||||
% feedback response (ctrl_null) against the stabilizing operation-mode
|
||||
% controller (ctrl_operation) holding T_avg at the nominal setpoint.
|
||||
%
|
||||
% Controller signature (thesis naming):
|
||||
% u = ctrl_<mode>(t, x, plant, ref)
|
||||
% where x is the continuous state vector and u is the external rod reactivity.
|
||||
|
||||
clear; clc; close all;
|
||||
addpath('controllers');
|
||||
|
||||
%% ===== Plant =====
|
||||
plant = pke_params();
|
||||
|
||||
%% ===== Scenario =====
|
||||
|
||||
% SG heat removal: step down to 80% at t = 30s, hold.
|
||||
Q_sg = @(t) plant.P0 * (1.0 - 0.2 * (t >= 30));
|
||||
|
||||
% Reference for the operation-mode controller.
|
||||
ref = struct();
|
||||
ref.T_avg = plant.T_c0; % hold T_avg at steady-state value
|
||||
|
||||
% Simulation window.
|
||||
tspan = [0, 600];
|
||||
|
||||
%% ===== Baseline: unforced feedback response =====
|
||||
fprintf('\n----- Run 1: ctrl_null (plant feedback only) -----\n');
|
||||
[t1, X1, U1] = pke_solver(plant, Q_sg, @ctrl_null, [], tspan);
|
||||
|
||||
%% ===== Closed loop: operation mode =====
|
||||
fprintf('\n----- Run 2: ctrl_operation (P on T_avg) -----\n');
|
||||
[t2, X2, U2] = pke_solver(plant, Q_sg, @ctrl_operation, ref, tspan);
|
||||
|
||||
%% ===== Plots =====
|
||||
plot_pke_results(t1, X1, U1, plant, Q_sg, 'ctrl_null (baseline)');
|
||||
plot_pke_results(t2, X2, U2, plant, Q_sg, 'ctrl_operation (P on T_avg)');
|
||||
|
||||
%% ===== Quick comparison figure =====
|
||||
CtoF = @(T) T*9/5 + 32;
|
||||
figure('Position', [100 50 1000 450], 'Name', 'T_avg comparison');
|
||||
plot(t1, CtoF(X1(:,9)), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t2, CtoF(X2(:,9)), 'r-', 'LineWidth', 1.5);
|
||||
yline(CtoF(plant.T_c0), 'k--', 'LineWidth', 1.0);
|
||||
xlabel('Time [s]'); ylabel('T_{avg} [F]');
|
||||
legend('ctrl\_null', 'ctrl\_operation', 'setpoint', 'Location', 'best');
|
||||
title('T_{avg} tracking: baseline vs stabilizing controller');
|
||||
grid on;
|
||||
@ -1,138 +0,0 @@
|
||||
%% main_mode_sweep.m — demo each DRC continuous-mode controller
|
||||
%
|
||||
% One scenario per mode, from an initial condition plausible for that
|
||||
% mode. Figures saved to ../docs/figures/mode_sweep_*.png so the runs
|
||||
% are reviewable without the MATLAB IDE.
|
||||
%
|
||||
% Modes, in the order DRC visits them starting from shutdown:
|
||||
% 1. shutdown — deep subcritical, cold IC, Q_sg = 0
|
||||
% 2. heatup — same cold IC, ramped T_avg reference
|
||||
% 3. operation — operating IC, 100% -> 80% Q_sg step
|
||||
% 4. scram — operating IC, rods slammed in, Q_sg decays to 0
|
||||
%
|
||||
% Each run produces plot_pke_results 4-panel plus saves a lightweight
|
||||
% summary figure.
|
||||
|
||||
clear; clc; close all;
|
||||
addpath('controllers');
|
||||
addpath(fullfile('..', 'reachability')); % for load_predicates
|
||||
|
||||
plant = pke_params();
|
||||
pred = load_predicates(plant); % T_standby and FRET-predicate concretizations
|
||||
figdir = fullfile('..', 'docs', 'figures');
|
||||
if ~exist(figdir, 'dir'), mkdir(figdir); end
|
||||
|
||||
CtoF = @(T) T*9/5 + 32;
|
||||
|
||||
%% ===== IC helpers =====
|
||||
% Operating IC: full-power steady state.
|
||||
x0_op = pke_initial_conditions(plant);
|
||||
|
||||
% Hot-standby shutdown IC: coolant flat at T_standby = T_c0 - 60 F ~ 275 C,
|
||||
% reactor deep subcritical with only a trace neutron population. T_standby
|
||||
% comes from reachability/predicates.json (single source of truth). Well
|
||||
% inside the model's +/-50 C trust region around the operating point, and
|
||||
% above coolant saturation at reduced operating pressure.
|
||||
T_standby = pred.constants.T_standby;
|
||||
n_shut = 1e-6;
|
||||
C_shut = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_shut;
|
||||
x0_shut = [n_shut; C_shut; T_standby; T_standby; T_standby];
|
||||
|
||||
% Heatup IC: reactor already taken critical at 0.1% power at hot-standby
|
||||
% temperature. Mirrors real plant practice: achieve criticality, then
|
||||
% heat up. Same T_standby as the shutdown IC — heatup begins from where
|
||||
% shutdown left off.
|
||||
n_heat = 1e-3;
|
||||
C_heat = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n_heat;
|
||||
x0_heat = [n_heat; C_heat; T_standby; T_standby; T_standby];
|
||||
|
||||
%% ===== Mode 1: SHUTDOWN =====
|
||||
fprintf('\n===== Mode 1: ctrl_shutdown =====\n');
|
||||
Q_sg_shut = @(t) 0; % no SG demand
|
||||
tspan_shut = [0, 600];
|
||||
[t1, X1, U1] = pke_solver(plant, Q_sg_shut, @ctrl_shutdown, [], tspan_shut, x0_shut);
|
||||
|
||||
%% ===== Mode 2: HEATUP =====
|
||||
fprintf('\n===== Mode 2: ctrl_heatup =====\n');
|
||||
ref_heatup = struct();
|
||||
ref_heatup.T_start = T_standby; % ~275 C (hot standby, = T_c0 - 60 F)
|
||||
ref_heatup.T_target = plant.T_c0; % 308.35 C (operating setpoint)
|
||||
ref_heatup.ramp_rate = 28 / 3600; % 28 C/hr, tech-spec limit
|
||||
Q_sg_heat = @(t) 0; % no SG demand during heatup
|
||||
tspan_heat = [0, 5400]; % ~90 min (33 C span at 28 C/hr = 71 min + settling)
|
||||
[t2, X2, U2] = pke_solver(plant, Q_sg_heat, @ctrl_heatup, ref_heatup, tspan_heat, x0_heat);
|
||||
|
||||
%% ===== Mode 3a: OPERATION (plain P) =====
|
||||
fprintf('\n===== Mode 3a: ctrl_operation (P) =====\n');
|
||||
Q_sg_op = @(t) plant.P0 * (1.0 - 0.2 * (t >= 30)); % 100% -> 80% step
|
||||
ref_op = struct('T_avg', plant.T_c0);
|
||||
tspan_op = [0, 600];
|
||||
[t3, X3, U3] = pke_solver(plant, Q_sg_op, @ctrl_operation, ref_op, tspan_op, x0_op);
|
||||
|
||||
%% ===== Mode 3b: OPERATION (LQR) =====
|
||||
fprintf('\n===== Mode 3b: ctrl_operation_lqr =====\n');
|
||||
[t3b, X3b, U3b] = pke_solver(plant, Q_sg_op, @ctrl_operation_lqr, [], tspan_op, x0_op);
|
||||
|
||||
%% ===== Mode 4: SCRAM =====
|
||||
fprintf('\n===== Mode 4: ctrl_scram =====\n');
|
||||
% After a scram signal, the turbine trips and SG isolation occurs fast.
|
||||
% Q_sg drops from P0 to a decay-heat-level sink (3% of P0) in ~10 s,
|
||||
% then holds. Not a true decay-heat model — good enough to show the
|
||||
% post-scram thermal trajectory without coolant going unphysically cold.
|
||||
Q_sg_scr = @(t) plant.P0 * max(0.03, 1 - max(0, t - 10) / 10);
|
||||
tspan_scr = [0, 600];
|
||||
[t4, X4, U4] = pke_solver(plant, Q_sg_scr, @ctrl_scram, [], tspan_scr, x0_op);
|
||||
|
||||
%% ===== Per-mode 4-panel plots =====
|
||||
plot_pke_results(t1, X1, U1, plant, Q_sg_shut, 'ctrl\_shutdown (cold IC)');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_1_shutdown.png'), 'Resolution', 150);
|
||||
|
||||
plot_pke_results(t2, X2, U2, plant, Q_sg_heat, 'ctrl\_heatup (ramp T\_avg)');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_2_heatup.png'), 'Resolution', 150);
|
||||
|
||||
plot_pke_results(t3, X3, U3, plant, Q_sg_op, 'ctrl\_operation (P on T\_avg)');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3_operation.png'), 'Resolution', 150);
|
||||
|
||||
plot_pke_results(t3b, X3b, U3b, plant, Q_sg_op, 'ctrl\_operation\_lqr');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_3b_operation_lqr.png'), 'Resolution', 150);
|
||||
|
||||
plot_pke_results(t4, X4, U4, plant, Q_sg_scr, 'ctrl\_scram');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_4_scram.png'), 'Resolution', 150);
|
||||
|
||||
%% ===== Heatup ramp-tracking figure =====
|
||||
% Overlay the T_ref signal on T_avg so the lag is visible.
|
||||
T_ref_trace = min(ref_heatup.T_start + ref_heatup.ramp_rate .* t2, ref_heatup.T_target);
|
||||
figure('Position', [100 80 900 400], 'Name', 'Heatup ramp tracking');
|
||||
plot(t2/60, CtoF(T_ref_trace), 'k--', 'LineWidth', 1.2); hold on;
|
||||
plot(t2/60, CtoF(X2(:,9)), 'r-', 'LineWidth', 1.5);
|
||||
xlabel('Time [min]'); ylabel('T_{avg} [F]');
|
||||
legend('T_{ref}(t) (28 C/hr ramp)', 'T_{avg} (plant)', 'Location', 'southeast');
|
||||
title('ctrl\_heatup: reference tracking');
|
||||
grid on;
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_heatup_tracking.png'), 'Resolution', 150);
|
||||
|
||||
%% ===== P vs LQR head-to-head on T_avg =====
|
||||
figure('Position', [100 80 900 400], 'Name', 'Operation: P vs LQR');
|
||||
plot(t3, CtoF(X3(:,9)), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t3b, CtoF(X3b(:,9)), 'r-', 'LineWidth', 1.5);
|
||||
yline(CtoF(plant.T_c0), 'k--');
|
||||
xlabel('Time [s]'); ylabel('T_{avg} [F]');
|
||||
legend('ctrl\_operation (P)', 'ctrl\_operation\_lqr', 'setpoint', 'Location', 'best');
|
||||
title('Operation mode: P vs LQR under 100% \rightarrow 80% Q_{sg} step');
|
||||
grid on;
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_op_P_vs_LQR.png'), 'Resolution', 150);
|
||||
|
||||
%% ===== Summary figure: normalized power across all modes =====
|
||||
figure('Position', [100 80 1000 450], 'Name', 'Normalized power, all modes');
|
||||
subplot(2,2,1); semilogy(t1, max(X1(:,1), 1e-20)); grid on;
|
||||
title('Shutdown'); xlabel('t [s]'); ylabel('n (log)');
|
||||
subplot(2,2,2); plot(t2/60, X2(:,1)); grid on;
|
||||
title('Heatup'); xlabel('t [min]'); ylabel('n');
|
||||
subplot(2,2,3); plot(t3, X3(:,1)); grid on;
|
||||
title('Operation'); xlabel('t [s]'); ylabel('n');
|
||||
subplot(2,2,4); semilogy(t4, max(X4(:,1), 1e-20)); grid on;
|
||||
title('Scram'); xlabel('t [s]'); ylabel('n (log)');
|
||||
sgtitle('Normalized reactor power n(t) across DRC modes');
|
||||
exportgraphics(gcf, fullfile(figdir, 'mode_sweep_power_overview.png'), 'Resolution', 150);
|
||||
|
||||
fprintf('\n=== Figures written to %s ===\n', figdir);
|
||||
@ -1,12 +0,0 @@
|
||||
function x0 = pke_initial_conditions(plant)
|
||||
% PKE_INITIAL_CONDITIONS Steady-state initial condition vector.
|
||||
%
|
||||
% Returns x0 = [n0; C1..C6; T_f0; T_c0; T_cold0] (10 x 1)
|
||||
% Precursor concentrations are the analytic equilibrium values at n = 1.
|
||||
|
||||
n0 = 1.0;
|
||||
C0 = (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * n0;
|
||||
|
||||
x0 = [n0; C0; plant.T_f0; plant.T_c0; plant.T_cold0];
|
||||
|
||||
end
|
||||
@ -1,54 +0,0 @@
|
||||
function [A, B, B_w, x_star, u_star, Q_star] = pke_linearize(plant, x_star, u_star, Q_star)
|
||||
% PKE_LINEARIZE Numerical Jacobians of the PKE+T/H RHS at an operating point.
|
||||
%
|
||||
% Computes (A, B, B_w) such that, for small perturbations dx, du, dw
|
||||
% around (x_star, u_star, Q_star):
|
||||
%
|
||||
% dx/dt ~ A*dx + B*du + B_w*dw, where w = Q_sg.
|
||||
%
|
||||
% Uses central finite differences. Step sizes are scaled to each
|
||||
% variable's magnitude; fall back to an absolute epsilon near zero.
|
||||
%
|
||||
% Inputs:
|
||||
% plant - parameter struct from pke_params()
|
||||
% x_star - (optional) 10x1 linearization point; default = operating x0
|
||||
% u_star - (optional) scalar control; default = 0
|
||||
% Q_star - (optional) scalar SG heat removal [W]; default = plant.P0
|
||||
%
|
||||
% Outputs:
|
||||
% A - 10x10 state Jacobian df/dx
|
||||
% B - 10x1 input Jacobian df/du
|
||||
% B_w - 10x1 disturbance Jacobian df/d(Q_sg)
|
||||
% x_star, u_star, Q_star - echoed back for the caller's records
|
||||
|
||||
if nargin < 2 || isempty(x_star), x_star = pke_initial_conditions(plant); end
|
||||
if nargin < 3 || isempty(u_star), u_star = 0; end
|
||||
if nargin < 4 || isempty(Q_star), Q_star = plant.P0; end
|
||||
|
||||
n = numel(x_star);
|
||||
eps_rel = 1e-6;
|
||||
eps_abs = 1e-8;
|
||||
|
||||
% Wrap the RHS as a function of (x, u, w) with plant fixed.
|
||||
f = @(x, u, w) pke_th_rhs(0, x, plant, @(t) w, u);
|
||||
|
||||
f0 = f(x_star, u_star, Q_star); %#ok<NASGU> % sanity slot
|
||||
|
||||
% --- A: df/dx via central differences, column by column ---
|
||||
A = zeros(n, n);
|
||||
for k = 1:n
|
||||
h = max(eps_rel * abs(x_star(k)), eps_abs);
|
||||
xp = x_star; xp(k) = xp(k) + h;
|
||||
xm = x_star; xm(k) = xm(k) - h;
|
||||
A(:, k) = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) / (2*h);
|
||||
end
|
||||
|
||||
% --- B: df/du ---
|
||||
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)) / (2*h);
|
||||
|
||||
% --- B_w: df/dQ_sg ---
|
||||
h = max(eps_rel * abs(Q_star), 1.0); % Q_star ~ 1e9, abs floor at 1 W
|
||||
B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) / (2*h);
|
||||
|
||||
end
|
||||
@ -1,86 +0,0 @@
|
||||
function p = pke_params()
|
||||
% PKE_PARAMS Returns the parameter struct for the coupled PKE + T/H model.
|
||||
%
|
||||
% All plant parameters, kinetic data, and steady-state conditions live here.
|
||||
% Modify this file to change the reactor you're modeling.
|
||||
|
||||
p = struct();
|
||||
|
||||
% --- Neutronics ---
|
||||
% Prompt neutron generation time for a typical PWR with enriched UO2 fuel
|
||||
% and a water-moderated thermal spectrum. Literature range: 1e-5 to 5e-4 s;
|
||||
% 1e-4 is a standard mid-range value (Duderstadt & Hamilton, Ch. 6).
|
||||
p.Lambda = 1e-4; % prompt neutron generation time [s]
|
||||
|
||||
% 6-group delayed neutron data for U-235 thermal fission (Keepin et al.,
|
||||
% Phys. Rev. 107, 1957). These are the canonical values used in virtually
|
||||
% all PKE models. Total beta_eff ~ 650 pcm.
|
||||
p.beta_i = [0.000215; 0.001424; 0.001274; 0.002568; 0.000748; 0.000273];
|
||||
p.lambda_i = [0.0124; 0.0305; 0.111; 0.301; 1.14; 3.01];
|
||||
p.beta = sum(p.beta_i);
|
||||
|
||||
% --- Thermal-hydraulic ---
|
||||
% 1000 MWth is a round number representative of a 2-loop PWR (e.g.,
|
||||
% Westinghouse 2-loop class). Actual plants range ~1800-3400 MWth for
|
||||
% 2-loop to 4-loop designs.
|
||||
p.P0 = 1000e6; % nominal thermal power [W]
|
||||
|
||||
% Fuel mass: representative of total UO2 inventory in the core.
|
||||
% A typical 4-loop PWR has ~100,000 kg UO2; scaled down for 1000 MWth.
|
||||
p.M_f = 50000; % fuel mass [kg]
|
||||
|
||||
% UO2 specific heat at operating temperatures (~400-800 C). Published
|
||||
% values range 250-350 J/kg-K (MATPRO / Fink correlation); 300 is a
|
||||
% standard lumped value for transient analysis.
|
||||
p.c_f = 300; % fuel specific heat [J/kg-K]
|
||||
|
||||
% Coolant mass in the core region. Based on core volume fraction of
|
||||
% water in a PWR (~40% moderator-to-fuel ratio) and active core volume.
|
||||
p.M_c = 20000; % coolant mass in core [kg]
|
||||
|
||||
% Pressurized water specific heat at ~15.5 MPa, ~300 C. NIST/IAPWS
|
||||
% gives cp ~ 5.4-5.5 kJ/kg-K in this range.
|
||||
p.c_c = 5450; % coolant specific heat [J/kg-K]
|
||||
|
||||
% Lumped fuel-to-coolant heat transfer coefficient x area. Chosen so that
|
||||
% the steady-state fuel-to-coolant dT = P0/hA = 20 C, which is in the
|
||||
% right ballpark for the average pellet-to-bulk-coolant temperature drop
|
||||
% in a lumped single-node fuel model.
|
||||
p.hA = 5e7; % fuel-to-coolant heat transfer coeff * area [W/K]
|
||||
|
||||
% Primary coolant mass flow rate. Typical PWR: ~4000-6000 kg/s per loop.
|
||||
% 5000 kg/s gives core dT = P0/(W*c_c) ~ 36.7 C, consistent with a
|
||||
% ~35-38 C core rise in operating PWRs.
|
||||
p.W = 5000; % coolant mass flow rate [kg/s]
|
||||
|
||||
% Coolant mass in SG + hot/cold leg piping. Represents the ex-core primary
|
||||
% inventory that participates in the loop transport delay. Larger than M_c
|
||||
% because the SG tubes and piping hold significant volume.
|
||||
p.M_sg = 30000; % coolant mass in SG + hot/cold leg piping [kg]
|
||||
|
||||
% --- Reactivity feedback coefficients ---
|
||||
% Doppler coefficient: fuel temperature feedback, predominantly from
|
||||
% U-238 resonance broadening. Typical BOL range: -2 to -3 pcm/K
|
||||
% (Todreas & Kazimi). -2.5 pcm/K is a standard mid-cycle value.
|
||||
p.alpha_f = -2.5e-5; % fuel (Doppler) [dk/k per K]
|
||||
|
||||
% Moderator temperature coefficient: captures water density change with
|
||||
% temperature and its effect on neutron moderation/absorption. Typical
|
||||
% range: -5 to -50 pcm/K depending on boron concentration and burnup.
|
||||
% -10 pcm/K is moderate (mid-cycle, moderate boron).
|
||||
p.alpha_c = -1.0e-4; % coolant (moderator) [dk/k per K]
|
||||
|
||||
% --- Steady-state (derived from energy balance at P0, n=1) ---
|
||||
% Q_sg = P0, all dT/dt = 0
|
||||
% T_hot - T_cold = P0 / (W*c_c)
|
||||
% T_c = (T_hot + T_cold) / 2
|
||||
% T_f = T_c + P0 / hA
|
||||
% Cold leg temperature of 290 C (554 F) is typical for a PWR at full
|
||||
% power (Westinghouse plants: ~288-293 C depending on design).
|
||||
p.T_cold0 = 290; % inlet coolant [C]
|
||||
p.dT_core = p.P0 / (p.W * p.c_c); % core delta-T [C]
|
||||
p.T_hot0 = p.T_cold0 + p.dT_core;
|
||||
p.T_c0 = (p.T_hot0 + p.T_cold0) / 2; % avg coolant [C]
|
||||
p.T_f0 = p.T_c0 + p.P0 / p.hA; % fuel [C]
|
||||
|
||||
end
|
||||
@ -1,64 +0,0 @@
|
||||
function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan, x0)
|
||||
% PKE_SOLVER Solve the coupled PKE + T/H system in closed loop with a mode controller.
|
||||
%
|
||||
% Inputs:
|
||||
% plant - parameter struct from pke_params()
|
||||
% Q_sg - function handle Q_sg(t) returning SG heat removal [W]
|
||||
% ctrl_fn - function handle u = ctrl_fn(t, x, plant, ref)
|
||||
% ref - struct of per-mode setpoints (e.g. ref.T_avg); [] if unused
|
||||
% tspan - [t_start, t_end] in seconds
|
||||
% x0 - (optional) initial state vector (10 x 1). Defaults to the
|
||||
% operating steady state from pke_initial_conditions(plant).
|
||||
%
|
||||
% Outputs:
|
||||
% t - time vector (M x 1) from the ODE solver
|
||||
% X - state matrix (M x 10): columns are [n, C1..C6, T_f, T_c, T_cold]
|
||||
% U - control input trajectory (M x 1): u evaluated at each (t(k), X(k,:))
|
||||
|
||||
CtoF = @(T) T * 9/5 + 32;
|
||||
|
||||
%% ===== Print Steady-State =====
|
||||
fprintf('=== Steady-State Conditions ===\n');
|
||||
fprintf(' beta = %.5f (%.1f pcm)\n', plant.beta, plant.beta*1e5);
|
||||
fprintf(' Lambda = %.1e s\n', plant.Lambda);
|
||||
fprintf(' P0 = %.0f MWth\n', plant.P0/1e6);
|
||||
fprintf(' T_cold = %.1f F\n', CtoF(plant.T_cold0));
|
||||
fprintf(' T_avg = %.1f F\n', CtoF(plant.T_c0));
|
||||
fprintf(' T_hot = %.1f F\n', CtoF(plant.T_hot0));
|
||||
fprintf(' T_fuel = %.1f F\n', CtoF(plant.T_f0));
|
||||
fprintf(' Core dT = %.1f F\n', plant.dT_core * 9/5);
|
||||
fprintf('================================\n\n');
|
||||
|
||||
%% ===== Solve closed-loop =====
|
||||
if nargin < 6 || isempty(x0)
|
||||
x0 = pke_initial_conditions(plant);
|
||||
end
|
||||
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
|
||||
|
||||
rhs = @(t, x) pke_th_rhs(t, x, plant, Q_sg, ctrl_fn(t, x, plant, ref));
|
||||
[t, X] = ode15s(rhs, tspan, x0, options);
|
||||
|
||||
%% ===== Reconstruct control trajectory =====
|
||||
M = length(t);
|
||||
U = zeros(M, 1);
|
||||
for k = 1:M
|
||||
U(k) = ctrl_fn(t(k), X(k,:).', plant, ref);
|
||||
end
|
||||
|
||||
%% ===== Print Final State =====
|
||||
n = X(end, 1); T_f = X(end, 8); T_c = X(end, 9); T_cold = X(end, 10);
|
||||
T_hot = 2*T_c - T_cold;
|
||||
rho_tot = U(end) ...
|
||||
+ plant.alpha_f*(T_f - plant.T_f0) ...
|
||||
+ plant.alpha_c*(T_c - plant.T_c0);
|
||||
|
||||
fprintf('=== Final State (t = %.0f s) ===\n', t(end));
|
||||
fprintf(' Power = %.1f MWth (%.3f x nominal)\n', n*plant.P0/1e6, n);
|
||||
fprintf(' T_fuel = %.1f F\n', CtoF(T_f));
|
||||
fprintf(' T_hot = %.1f F\n', CtoF(T_hot));
|
||||
fprintf(' T_avg = %.1f F\n', CtoF(T_c));
|
||||
fprintf(' T_cold = %.1f F\n', CtoF(T_cold));
|
||||
fprintf(' u = %.4f $ (external reactivity from controller)\n', U(end)/plant.beta);
|
||||
fprintf(' rho = %.4f $ (total, incl. T-feedback)\n', rho_tot/plant.beta);
|
||||
|
||||
end
|
||||
@ -1,49 +0,0 @@
|
||||
function dxdt = pke_th_rhs(t, x, plant, Q_sg, u)
|
||||
% PKE_TH_RHS ODE right-hand side for coupled PKE + thermal-hydraulics.
|
||||
%
|
||||
% State x = [n; C1..C6; T_f; T_c; T_cold] (10 x 1)
|
||||
%
|
||||
% Inputs:
|
||||
% t - time [s]
|
||||
% x - continuous state vector (10 x 1)
|
||||
% plant - parameter struct from pke_params()
|
||||
% Q_sg - function handle Q_sg(t) returning SG heat removal [W]
|
||||
% u - external reactivity (rod worth) [dk/k], scalar
|
||||
%
|
||||
% The controller is responsible for computing u; this function treats u
|
||||
% as a given input. Q_sg is the bounded disturbance.
|
||||
|
||||
n = x(1);
|
||||
C = x(2:7);
|
||||
T_f = x(8);
|
||||
T_c = x(9);
|
||||
T_cold = x(10);
|
||||
T_hot = 2 * T_c - T_cold;
|
||||
|
||||
Q_sg_val = Q_sg(t);
|
||||
|
||||
% --- Reactivity with feedback ---
|
||||
% rho = u + alpha_f*(T_f - T_f0) + alpha_c*(T_c - T_c0)
|
||||
rho = u ...
|
||||
+ plant.alpha_f * (T_f - plant.T_f0) ...
|
||||
+ plant.alpha_c * (T_c - plant.T_c0);
|
||||
|
||||
% --- Neutronics ---
|
||||
dndt = (rho - plant.beta) / plant.Lambda * n + plant.lambda_i' * C;
|
||||
dCdt = (plant.beta_i / plant.Lambda) * n - plant.lambda_i .* C;
|
||||
|
||||
% --- Thermal-hydraulics ---
|
||||
% Fuel node: fission power in, heat transfer to coolant out
|
||||
dTf_dt = (plant.P0 * n - plant.hA * (T_f - T_c)) / (plant.M_f * plant.c_f);
|
||||
|
||||
% Core coolant: heat from fuel in, bulk enthalpy flow out
|
||||
dTc_dt = (plant.hA * (T_f - T_c) - 2*plant.W*plant.c_c*(T_c - T_cold)) ...
|
||||
/ (plant.M_c * plant.c_c);
|
||||
|
||||
% SG/loop node: hot leg enthalpy in, SG heat removal out
|
||||
dTcold_dt = (plant.W * plant.c_c * (T_hot - T_cold) - Q_sg_val) ...
|
||||
/ (plant.M_sg * plant.c_c);
|
||||
|
||||
dxdt = [dndt; dCdt; dTf_dt; dTc_dt; dTcold_dt];
|
||||
|
||||
end
|
||||
@ -1,76 +0,0 @@
|
||||
function plot_pke_results(t, X, U, plant, Q_sg, plot_title)
|
||||
% PLOT_PKE_RESULTS Standard 4-panel plot of PKE + T/H simulation results.
|
||||
%
|
||||
% Inputs:
|
||||
% t - time vector from ODE solver
|
||||
% X - state matrix (M x 10): [n, C1..C6, T_f, T_c, T_cold]
|
||||
% U - control input trajectory (M x 1): scalar u at each time
|
||||
% plant - parameter struct from pke_params()
|
||||
% Q_sg - function handle Q_sg(t) used in the simulation
|
||||
% plot_title - (optional) string used in the figure name and first subplot
|
||||
|
||||
if nargin < 6 || isempty(plot_title)
|
||||
plot_title = 'PKE + T/H closed-loop';
|
||||
end
|
||||
|
||||
% --- Extract states ---
|
||||
n = X(:, 1);
|
||||
T_f = X(:, 8);
|
||||
T_c = X(:, 9);
|
||||
T_cold = X(:, 10);
|
||||
T_hot = 2 * T_c - T_cold;
|
||||
|
||||
% --- Convert temperatures to Fahrenheit for display ---
|
||||
CtoF = @(T) T * 9/5 + 32;
|
||||
T_f_F = CtoF(T_f);
|
||||
T_c_F = CtoF(T_c);
|
||||
T_cold_F = CtoF(T_cold);
|
||||
T_hot_F = CtoF(T_hot);
|
||||
|
||||
% --- Derived quantities ---
|
||||
power_MW = n * plant.P0 / 1e6;
|
||||
Q_sg_MW = arrayfun(Q_sg, t) / 1e6;
|
||||
|
||||
rho_fb = plant.alpha_f * (T_f - plant.T_f0) + plant.alpha_c * (T_c - plant.T_c0);
|
||||
rho_tot = U + rho_fb;
|
||||
|
||||
% --- Plot ---
|
||||
figure('Position', [100 50 1000 950], 'Name', plot_title);
|
||||
|
||||
% Power & SG demand
|
||||
subplot(4, 1, 1);
|
||||
plot(t, power_MW, 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t, Q_sg_MW, 'r--', 'LineWidth', 1.2);
|
||||
xlabel('Time [s]'); ylabel('Power [MW_{th}]');
|
||||
legend('Reactor Power', 'Q_{sg}', 'Location', 'east');
|
||||
title(sprintf('%s — Power and SG Heat Demand', plot_title));
|
||||
grid on;
|
||||
|
||||
% Reactivity: controller u, feedback, total
|
||||
subplot(4, 1, 2);
|
||||
plot(t, U / plant.beta, 'k-', 'LineWidth', 1.5); hold on;
|
||||
plot(t, rho_fb / plant.beta, 'b-.', 'LineWidth', 1.2);
|
||||
plot(t, rho_tot / plant.beta, 'm:', 'LineWidth', 1.5);
|
||||
xlabel('Time [s]'); ylabel('\rho [$]');
|
||||
legend('u (controller)', 'Feedback', 'Total', 'Location', 'east');
|
||||
title('Reactivity');
|
||||
grid on;
|
||||
|
||||
% Coolant temperatures
|
||||
subplot(4, 1, 3);
|
||||
plot(t, T_hot_F, 'r-', 'LineWidth', 1.5); hold on;
|
||||
plot(t, T_c_F, 'm--', 'LineWidth', 1.2);
|
||||
plot(t, T_cold_F, 'b-', 'LineWidth', 1.5);
|
||||
xlabel('Time [s]'); ylabel('Temperature [F]');
|
||||
legend('T_{hot}', 'T_{avg}', 'T_{cold}', 'Location', 'east');
|
||||
title('Coolant Temperatures');
|
||||
grid on;
|
||||
|
||||
% Fuel temperature
|
||||
subplot(4, 1, 4);
|
||||
plot(t, T_f_F, 'Color', [0.85 0.33 0.1], 'LineWidth', 1.5);
|
||||
xlabel('Time [s]'); ylabel('Temperature [F]');
|
||||
title('Fuel Temperature');
|
||||
grid on;
|
||||
|
||||
end
|
||||
@ -1,20 +0,0 @@
|
||||
%% test_heatup_rate.m — measure the actual heatup-rate trajectory from ctrl_heatup
|
||||
clear; addpath('controllers'); addpath('../reachability');
|
||||
plant = pke_params();
|
||||
pred = load_predicates(plant);
|
||||
x0 = [1e-3; (plant.beta_i ./ (plant.lambda_i * plant.Lambda)) * 1e-3; ...
|
||||
pred.constants.T_standby; pred.constants.T_standby; pred.constants.T_standby];
|
||||
ref = struct('T_start', pred.constants.T_standby, ...
|
||||
'T_target', plant.T_c0, ...
|
||||
'ramp_rate', 28/3600);
|
||||
[t, X, U] = pke_solver(plant, @(t) 0, @ctrl_heatup, ref, [0 5400], x0);
|
||||
|
||||
a_f = 0.4587; a_c = -0.9587; a_cold = 0.5000;
|
||||
dTcdt = a_f*X(:,8) + a_c*X(:,9) + a_cold*X(:,10); % C/s
|
||||
|
||||
fprintf('\n=== Heatup rate statistics ===\n');
|
||||
fprintf(' max dT_c/dt: %.4f C/s = %6.1f C/hr\n', max(dTcdt), max(dTcdt)*3600);
|
||||
fprintf(' min dT_c/dt: %.4f C/s = %6.1f C/hr\n', min(dTcdt), min(dTcdt)*3600);
|
||||
fprintf(' rate limit: +/- 0.01389 C/s = +/- 50.0 C/hr\n');
|
||||
fprintf(' violates upper? %s\n', string(any(dTcdt > 0.01389)));
|
||||
fprintf(' violates lower? %s\n', string(any(dTcdt < -0.01389)));
|
||||
@ -1,106 +0,0 @@
|
||||
%% test_linearize.m — sanity-check the Jacobians against the nonlinear sim
|
||||
%
|
||||
% Simulate two parallel trajectories under a small Q_sg step (5% down):
|
||||
% 1. Full nonlinear plant with u = 0
|
||||
% 2. Linear dx/dt = A*dx + B_w*dw propagated by ode45
|
||||
% Compare deviations from x_op. For a small-enough disturbance, the
|
||||
% trajectories should overlay to within a few percent over a few minutes.
|
||||
|
||||
clear; clc; close all;
|
||||
addpath('controllers');
|
||||
|
||||
plant = pke_params();
|
||||
x_op = pke_initial_conditions(plant);
|
||||
|
||||
[A, B, B_w, xs, us, Qs] = pke_linearize(plant, x_op, 0, plant.P0);
|
||||
|
||||
%% ===== Eigenvalue/conditioning check =====
|
||||
eigs_A = eig(A);
|
||||
fprintf('\n=== Linearization at operating point ===\n');
|
||||
fprintf(' ||A|| = %.3e\n', norm(A));
|
||||
fprintf(' ||B|| = %.3e\n', norm(B));
|
||||
fprintf(' ||B_w|| = %.3e\n', norm(B_w));
|
||||
fprintf(' stable? = %s (max Re(eig) = %.3e)\n', ...
|
||||
iif(all(real(eigs_A) < 1e-8), 'YES', 'NO'), max(real(eigs_A)));
|
||||
fprintf(' fastest pole ~ %.3e /s\n', max(abs(real(eigs_A))));
|
||||
fprintf(' slowest pole ~ %.3e /s\n', min(abs(real(eigs_A(real(eigs_A) < -1e-10)))));
|
||||
|
||||
%% ===== Linear vs nonlinear trajectory under small Q_sg step =====
|
||||
% 5% down step at t = 30s.
|
||||
dQ = -0.05 * plant.P0;
|
||||
Q_sg = @(t) plant.P0 + dQ * (t >= 30);
|
||||
|
||||
tspan = [0, 300];
|
||||
|
||||
% Nonlinear
|
||||
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
|
||||
rhs_nl = @(t, x) pke_th_rhs(t, x, plant, Q_sg, 0);
|
||||
[t_nl, X_nl] = ode15s(rhs_nl, tspan, x_op, opts);
|
||||
|
||||
% Linear (deviation)
|
||||
dQ_of_t = @(t) dQ * (t >= 30);
|
||||
rhs_lin = @(t, dx) A*dx + B_w*dQ_of_t(t);
|
||||
[t_lin, DX_lin] = ode15s(rhs_lin, tspan, zeros(size(x_op)), opts);
|
||||
|
||||
% Reconstruct absolute linear state for plotting
|
||||
X_lin = DX_lin + x_op.';
|
||||
|
||||
%% ===== Plot comparison — focus on thermal states =====
|
||||
CtoF = @(T) T*9/5 + 32;
|
||||
figdir = fullfile('..', 'docs', 'figures');
|
||||
if ~exist(figdir, 'dir'), mkdir(figdir); end
|
||||
|
||||
figure('Position', [100 80 1100 700], 'Name', 'Linear vs nonlinear sanity');
|
||||
|
||||
subplot(2,2,1);
|
||||
plot(t_nl, X_nl(:,1), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t_lin, X_lin(:,1), 'r--','LineWidth', 1.5);
|
||||
xlabel('t [s]'); ylabel('n'); grid on;
|
||||
title('Normalized power'); legend('nonlinear', 'linear', 'Location', 'best');
|
||||
|
||||
subplot(2,2,2);
|
||||
plot(t_nl, CtoF(X_nl(:,8)), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t_lin, CtoF(X_lin(:,8)), 'r--','LineWidth', 1.5);
|
||||
xlabel('t [s]'); ylabel('T_f [F]'); grid on;
|
||||
title('Fuel temperature');
|
||||
|
||||
subplot(2,2,3);
|
||||
plot(t_nl, CtoF(X_nl(:,9)), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t_lin, CtoF(X_lin(:,9)), 'r--','LineWidth', 1.5);
|
||||
xlabel('t [s]'); ylabel('T_c [F]'); grid on;
|
||||
title('Avg coolant');
|
||||
|
||||
subplot(2,2,4);
|
||||
plot(t_nl, CtoF(X_nl(:,10)), 'b-', 'LineWidth', 1.5); hold on;
|
||||
plot(t_lin, CtoF(X_lin(:,10)),'r--','LineWidth', 1.5);
|
||||
xlabel('t [s]'); ylabel('T_{cold} [F]'); grid on;
|
||||
title('Cold-leg coolant');
|
||||
|
||||
sgtitle('Linear-model sanity check: 5% Q_{sg} down-step, u = 0');
|
||||
exportgraphics(gcf, fullfile(figdir, 'linearize_sanity.png'), 'Resolution', 150);
|
||||
|
||||
%% ===== Quantitative error at a few times =====
|
||||
fprintf('\n=== Max |linear - nonlinear| at sampled times ===\n');
|
||||
fprintf(' (normalized by operating-point magnitude where nonzero)\n');
|
||||
for ts = [60, 120, 300]
|
||||
xi_nl = interp1(t_nl, X_nl, ts);
|
||||
xi_lin = interp1(t_lin, X_lin, ts);
|
||||
rel_err = abs(xi_nl - xi_lin) ./ max(abs(x_op.'), 1e-6);
|
||||
fprintf(' t = %3d s: max rel err = %.2e (%s)\n', ts, max(rel_err), ...
|
||||
state_name(find(rel_err == max(rel_err), 1)));
|
||||
end
|
||||
|
||||
%% ===== Save linearization =====
|
||||
save(fullfile('..', 'reachability', 'linearization_at_op.mat'), ...
|
||||
'A', 'B', 'B_w', 'xs', 'us', 'Qs', '-v7');
|
||||
fprintf('\nSaved A/B/B_w to ../reachability/linearization_at_op.mat\n');
|
||||
|
||||
%% ===== helpers =====
|
||||
function s = state_name(k)
|
||||
names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'};
|
||||
s = names{k};
|
||||
end
|
||||
|
||||
function y = iif(cond, a, b)
|
||||
if cond, y = a; else, y = b; end
|
||||
end
|
||||
@ -4,95 +4,95 @@ Continuous-mode verification for the PWR_HYBRID_3 hybrid controller.
|
||||
|
||||
## Soundness status: APPROXIMATE
|
||||
|
||||
The current `reach_operation.m` result is **not a sound reach tube for
|
||||
the physical plant**. It is a sound over-approximation of the
|
||||
*linearized* closed-loop system (A_cl = A - BK around x_op) under
|
||||
bounded disturbance. The linear model is itself an approximation of
|
||||
the nonlinear plant (`../plant-model/pke_th_rhs.m`), and that
|
||||
approximation error is not currently bounded or inflated into the tube.
|
||||
The current linear-reach results are **not a sound reach tube for
|
||||
the physical plant**. They are sound over-approximations of the
|
||||
*linearized* closed-loop system ($A_{\mathrm{cl}} = A - BK$ around
|
||||
$x_{\mathrm{op}}$) under bounded disturbance. The linear model is
|
||||
itself an approximation of the nonlinear plant (`../code/src/pke_th_rhs.jl`),
|
||||
and that approximation error is not currently bounded or inflated
|
||||
into the tube.
|
||||
|
||||
Two paths to upgrade to a sound result:
|
||||
|
||||
1. **Nonlinear reach directly** — CORA `nonlinearSys`, JuliaReach
|
||||
`BlackBoxContinuousSystem`, or equivalent. More expensive but the
|
||||
honest answer.
|
||||
1. **Nonlinear reach directly** — TMJets (Taylor-model integration)
|
||||
via ReachabilityAnalysis.jl. Currently limited to ~10-second
|
||||
horizons by prompt-neutron stiffness; needs a reduced-order PKE
|
||||
(prompt-jump approximation) to extend to mode-obligation horizons.
|
||||
2. **Linear reach + Taylor-remainder inflation** — compute an upper
|
||||
bound on `||f_nl(x, u) - (A x + B u)||` over the reach set (via
|
||||
Hessian norm estimate on each component of `f_nl`) and inflate the
|
||||
linear tube by that bound. Less expensive, still rigorous.
|
||||
bound on `||f_nl(x, u) - (A x + B u)||` over the reach set and
|
||||
inflate the linear tube by that bound. Cheaper, still rigorous.
|
||||
|
||||
Both are thesis-blocking for any safety claim. Deferred only until
|
||||
the per-mode plumbing is solid; it is not a "nice to have".
|
||||
|
||||
The current 5-orders-of-margin buffer (reach envelope ~0.03 K against
|
||||
a 5 K safety band) means linearization error would have to be huge to
|
||||
invalidate the conclusion, but that is vibes, not a proof.
|
||||
Both are thesis-blocking for any safety claim. The current
|
||||
5-orders-of-margin buffer (reach envelope ~0.03 K against a
|
||||
5 K safety band on $T_c$) means linearization error would have to
|
||||
be huge to invalidate the conclusion — but that's vibes, not a proof.
|
||||
|
||||
## Related open issues
|
||||
|
||||
- **Saturation semantics.** `ctrl_heatup.m` uses `sat(u, u_min, u_max)`.
|
||||
Saturation is formally a 3-mode piecewise-affine system. For
|
||||
heatup reach this has to be handled as (a) hybrid locations, or
|
||||
(b) proven dormant via reach on `u_unsat`. Not modeled in the
|
||||
- **Saturation semantics.** `ctrl_heatup.jl` uses `clamp(u, u_min, u_max)`.
|
||||
Saturation is formally a 3-mode piecewise-affine system. For
|
||||
heatup reach this must be handled as (a) hybrid locations, or
|
||||
(b) proven dormant via reach on `u_unsat`. Not modeled in the
|
||||
current artifacts (operation-mode LQR saturation is dormant in
|
||||
practice but the proof is implicit).
|
||||
- **Parametric uncertainty in α_f, α_c.** Real plants have α drift
|
||||
with burnup (~20%), boron (α_c ranges 10×), xenon. The
|
||||
feedback-linearization in `ctrl_heatup.m` assumes exact α; a robust
|
||||
- **Parametric uncertainty in α_f, α_c.** Real plants have α drift
|
||||
with burnup (~20%), boron (α_c ranges 10×), xenon. The
|
||||
feedback-linearization in `ctrl_heatup` assumes exact α; a robust
|
||||
treatment would make α an interval and propagate parametric reach.
|
||||
Currently idealized — flag in the chapter.
|
||||
Currently idealized.
|
||||
|
||||
## What's here
|
||||
|
||||
**Per-mode only.** Following the compositionality argument in the thesis:
|
||||
verify each continuous mode separately, let the DRC handle discrete
|
||||
switching. Current focus: **operation mode** under LQR feedback.
|
||||
**Per-mode only.** Following the compositionality argument in the
|
||||
thesis: verify each continuous mode separately, let the DRC handle
|
||||
discrete switching.
|
||||
|
||||
## What's here
|
||||
Files:
|
||||
|
||||
- `linearization_at_op.mat` — A, B, B_w and reference point, generated by
|
||||
`../plant-model/test_linearize.m`.
|
||||
- `reach_linear.m` — box-zonotope propagation of the closed-loop linear
|
||||
model under bounded disturbance. Pure MATLAB, no external toolbox.
|
||||
- `barrier_lyapunov.m` — Lyapunov-ellipsoid barrier certificate for the
|
||||
closed-loop linear system. Solves a Lyapunov equation, reports the
|
||||
smallest sub-level set containing the initial set and closed under
|
||||
the disturbance.
|
||||
- `reach_operation.m` — end-to-end operation-mode reach: linearize at
|
||||
x_op, compute LQR gain, propagate zonotope reach set, check against
|
||||
the `t_avg_in_range` predicate.
|
||||
- `figures/` — generated plots.
|
||||
- `predicates.json` — single source of truth for predicate
|
||||
concretizations. Three groups:
|
||||
- `operational_deadbands` — soft bands used by DRC for mode
|
||||
transitions (`t_avg_above_min`, `t_avg_in_range`, `p_above_crit`).
|
||||
- `safety_limits` — hard one-sided halfspaces (fuel centerline,
|
||||
trip setpoints, subcooling, heatup-rate bounds).
|
||||
- `mode_invariants` — `inv1_holds`, `inv2_holds` as conjunctions
|
||||
of safety limits.
|
||||
- `mode_boundaries` — per-mode $X_{\mathrm{entry}}$, $X_{\mathrm{safe}}$,
|
||||
$X_{\mathrm{exit}}$, $T_{\min}$, $T_{\max}$.
|
||||
|
||||
- `WALKTHROUGH.md` — standalone document explaining the reach-obligation
|
||||
taxonomy, per-mode entry/exit definitions, current results, soundness
|
||||
status. Read this for the full story.
|
||||
|
||||
The reach code itself lives in `../code/`:
|
||||
|
||||
- `../code/src/reach_linear.jl` — hand-rolled box reach propagator.
|
||||
- `../code/src/load_predicates.jl` — reads `predicates.json`.
|
||||
- `../code/scripts/reach_operation.jl` — operation-mode linear reach.
|
||||
- `../code/scripts/barrier_lyapunov.jl` — Lyapunov barrier attempt.
|
||||
- `../code/scripts/barrier_compare_OL_CL.jl` — OL vs CL comparison.
|
||||
- `../code/scripts/reach_heatup_nonlinear.jl` — TMJets nonlinear reach
|
||||
(heatup, saturation-disabled, 10-second horizon cap).
|
||||
|
||||
## Running
|
||||
|
||||
From MATLAB:
|
||||
|
||||
```matlab
|
||||
cd reachability
|
||||
reach_operation % computes reach set + plots
|
||||
barrier_lyapunov % solves Lyapunov, reports invariant ellipsoid
|
||||
```bash
|
||||
cd code
|
||||
julia --project=. scripts/reach_operation.jl
|
||||
julia --project=. scripts/barrier_lyapunov.jl
|
||||
julia --project=. scripts/barrier_compare_OL_CL.jl
|
||||
julia --project=. scripts/reach_heatup_nonlinear.jl
|
||||
```
|
||||
|
||||
## Tool choice
|
||||
|
||||
Currently using a hand-rolled zonotope reach because:
|
||||
- Avoids a ~0.5 GB CORA install for a first-pass result.
|
||||
- Linear reach with bounded disturbance has a clean analytic form
|
||||
(matrix exponential on the state, integral of e^(A(t-s))·B_w·w ds
|
||||
for the disturbance).
|
||||
- Stays inside MATLAB, which is where the plant model lives.
|
||||
|
||||
If we need nonlinear reach (and we will, for non-LQR controllers or
|
||||
larger reach sets where linearization error matters), the planned
|
||||
options are CORA (MATLAB) or JuliaReach (port the plant to Julia).
|
||||
Results save here as `*.mat` (gitignored).
|
||||
|
||||
## What this does NOT do yet
|
||||
|
||||
- Any sound reach tube (see top of this file).
|
||||
- Nonlinear reach for the original P controller on operation.
|
||||
- Heatup reach (ramped reference makes x* time-varying — needs
|
||||
trajectory-LQR or a different formulation, and the saturation
|
||||
semantics need to be made explicit).
|
||||
- Nonlinear reach on horizons > 10 s (needs reduced-order PKE).
|
||||
- Shutdown, scram, initialization reach.
|
||||
- Hybrid-system level verification (mode switching validity).
|
||||
- Parametric robustness to α_f, α_c drift.
|
||||
- Polytopic or SOS barriers — the canonical quadratic Lyapunov
|
||||
barrier fails structurally on this plant (see `WALKTHROUGH.md`
|
||||
and the OL-vs-CL comparison script).
|
||||
|
||||
@ -128,7 +128,13 @@ will require real tech-spec numbers.
|
||||
|
||||
## 3. Code walkthrough
|
||||
|
||||
Files in `reachability/` and what each does.
|
||||
Files in `code/` (Julia) and `reachability/` (JSON + docs) and what each does.
|
||||
|
||||
**Note:** as of 2026-04-20, all reach / plant / controller code is Julia.
|
||||
The pipeline was ported end-to-end from an earlier MATLAB implementation.
|
||||
Numerical results match MATLAB to machine precision on per-halfspace
|
||||
margins; trajectories match to ~3 decimals (different ODE solvers,
|
||||
same stiff family).
|
||||
|
||||
### `predicates.json`
|
||||
|
||||
@ -146,7 +152,7 @@ categories:
|
||||
|
||||
Plus `mode_boundaries` (entry/exit/time budgets per mode, as above).
|
||||
|
||||
### `load_predicates.m`
|
||||
### `code/src/load_predicates.jl`
|
||||
|
||||
Reads the JSON, resolves `rhs_expr` strings (which reference
|
||||
plant-derived constants like `T_c0`, `T_cold0`, `T_standby`) into numeric
|
||||
@ -163,7 +169,7 @@ pred.mode_invariants.inv2_holds % conjunction of safety_limits
|
||||
Each `A_poly` is an `n_halfspaces x 10` matrix; each `b_poly` is
|
||||
`n_halfspaces x 1`. The predicate is `{x : A_poly * x ≤ b_poly}`.
|
||||
|
||||
### `pke_linearize.m` (in `plant-model/`)
|
||||
### `code/src/pke_linearize.jl`
|
||||
|
||||
Numerical Jacobians via central finite differences at any
|
||||
`(x_star, u_star, Q_star)`. Returns `(A, B, B_w)` such that for small
|
||||
@ -173,7 +179,7 @@ perturbations `dx`, `du`, `dw`:
|
||||
dx/dt ≈ A dx + B du + B_w dw, w = Q_sg
|
||||
```
|
||||
|
||||
Validated in `plant-model/test_linearize.m`: under a 5% Q_sg step at
|
||||
Validated in `code/scripts/test_linearize.jl` (via `sim_sanity.jl`): under a 5% Q_sg step at
|
||||
`u=0`, the linear model matches the nonlinear sim to ~4e-4 relative
|
||||
error at 60 s, improving to 5e-6 at 300 s.
|
||||
|
||||
@ -189,7 +195,7 @@ See `docs/figures/linearize_sanity.png`:
|
||||
|
||||
|
||||
|
||||
### `reach_linear.m`
|
||||
### `code/src/reach_linear.jl`
|
||||
|
||||
~50 lines. Propagates a box over-approximation of the reach tube for a
|
||||
linear system `dx/dt = A_cl·x + B_w·w`, with `x(0)` in a box and `w` in
|
||||
@ -225,7 +231,7 @@ M_aug = expm([A_cl, B_w; zeros(1, n+1)] * dt);
|
||||
G_step = M_aug(1:n, n+1); % upper-right block IS the integral
|
||||
```
|
||||
|
||||
### `reach_operation.m`
|
||||
### `code/scripts/reach_operation.jl`
|
||||
|
||||
End-to-end operation-mode reach. Pipeline:
|
||||
|
||||
@ -255,9 +261,9 @@ for each halfspace a' x <= b in inv2_holds:
|
||||
|
||||
Reports per-halfspace margins so failure modes are attributable.
|
||||
|
||||
### `barrier_lyapunov.m`
|
||||
### `code/scripts/barrier_lyapunov.jl`
|
||||
|
||||
Analytic counterpart to `reach_operation.m`. Solves a Lyapunov equation
|
||||
Analytic counterpart to `code/scripts/reach_operation.jl`. Solves a Lyapunov equation
|
||||
to get an invariant ellipsoid, checks whether the ellipsoid fits inside
|
||||
the safety limits.
|
||||
|
||||
@ -277,10 +283,10 @@ max |a' dx| on {dx : dx' P dx ≤ γ} = sqrt(γ · a' P^(-1) a)
|
||||
```
|
||||
|
||||
Sweeps `Qbar(9,9)` from 10 to 1e6 to find the best quadratic barrier.
|
||||
See `barrier_compare_OL_CL.m` for the open-loop vs closed-loop
|
||||
See `code/scripts/barrier_compare_OL_CL.jl` for the open-loop vs closed-loop
|
||||
comparison.
|
||||
|
||||
### `barrier_compare_OL_CL.m`
|
||||
### `code/scripts/barrier_compare_OL_CL.jl`
|
||||
|
||||
Side-by-side: run the Lyapunov barrier on open-loop `A` (no controller)
|
||||
vs closed-loop `A_cl = A - BK`. Written to answer the question "is
|
||||
@ -361,14 +367,14 @@ it's a structural issue not a tuning one.
|
||||
|
||||
### Heatup mode — simulation only, no reach yet
|
||||
|
||||
The `main_mode_sweep.m` heatup scenario drives the plant from hot
|
||||
The `code/scripts/main_mode_sweep.jl` heatup scenario drives the plant from hot
|
||||
standby through the operating temperature:
|
||||
|
||||
|
||||
|
||||
Controller tracks the 28 °C/hr ramp with ~2 °F lag and ~5 °F final
|
||||
overshoot. **Peak rate during ramp is 48.5 °C/hr** (verified by
|
||||
`plant-model/test_heatup_rate.m`), right at the 50 °C/hr placeholder
|
||||
`code/scripts/test_heatup_rate.jl`), right at the 50 °C/hr placeholder
|
||||
invariant but 70% over the nominal tech-spec rate. A strict 28 °C/hr
|
||||
invariant would be violated by the current controller tuning.
|
||||
|
||||
@ -405,7 +411,7 @@ Power overview:
|
||||
|
||||
**The current reach tube is approximate, not sound, for the nonlinear plant.**
|
||||
|
||||
`reach_linear.m` produces a sound over-approximation of the *linearized*
|
||||
`code/src/reach_linear.jl` produces a sound over-approximation of the *linearized*
|
||||
closed-loop system's reach set. The linear model approximates the
|
||||
nonlinear plant with an error that is bounded but not currently
|
||||
computed or propagated into the tube. Two paths to soundness:
|
||||
@ -449,7 +455,7 @@ vibes, not a proof.
|
||||
### What's next
|
||||
|
||||
1. **Nonlinear reach on heatup via JuliaReach — partial result in.**
|
||||
`../julia-port/scripts/reach_heatup_nonlinear.jl` ran successfully
|
||||
`code/scripts/reach_heatup_nonlinear.jl` ran successfully
|
||||
to T=10 s on the full 10-state nonlinear closed-loop (saturation
|
||||
disabled). Got 10,583 reach-sets with envelope T_c ∈ [274.45, 295]
|
||||
and n ∈ [-5e-4, 5e-3]. Nonlinear reach is *functional* in Julia —
|
||||
|
||||
@ -1,86 +0,0 @@
|
||||
%% barrier_compare_OL_CL.m — does LQR help the Lyapunov barrier, or not?
|
||||
%
|
||||
% Head-to-head: solve the same Lyapunov-ellipsoid barrier analysis on
|
||||
% (i) the open-loop plant A (u = 0)
|
||||
% (ii) the LQR closed loop A_cl = A - B*K
|
||||
% with identical Qbar and identical disturbance bound. Report per-
|
||||
% halfspace bounds against inv2_holds so we can see exactly where LQR
|
||||
% helps and where it doesn't.
|
||||
%
|
||||
% The plant is already Hurwitz open-loop (max Re(eig) = -0.012 /s), so
|
||||
% the Lyapunov equation has a solution in both cases.
|
||||
|
||||
clear; clc;
|
||||
addpath('../plant-model', '../plant-model/controllers');
|
||||
|
||||
plant = pke_params();
|
||||
x_op = pke_initial_conditions(plant);
|
||||
pred = load_predicates(plant);
|
||||
|
||||
[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0);
|
||||
|
||||
%% ===== LQR gain =====
|
||||
Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1]);
|
||||
R_lqr = 1e6;
|
||||
try
|
||||
K = lqr(A, B, Q_lqr, R_lqr);
|
||||
catch
|
||||
[~, ~, K] = icare(A, B, Q_lqr, R_lqr);
|
||||
end
|
||||
A_cl = A - B*K;
|
||||
|
||||
%% ===== Shared pieces =====
|
||||
Qbar = diag([1, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1, 1e2, 1]); % best from the earlier sweep
|
||||
w_bar = 0.15 * plant.P0;
|
||||
|
||||
delta_entry = [0.01 * x_op(1);
|
||||
0.001 * abs(x_op(2:7));
|
||||
0.1; 0.1; 0.1];
|
||||
|
||||
inv2 = pred.mode_invariants.inv2_holds;
|
||||
A_inv = inv2.A_poly; b_inv = inv2.b_poly; comps = inv2.components;
|
||||
b_dev = b_inv - A_inv * x_op; % headroom in deviation coords
|
||||
|
||||
%% ===== Helper =====
|
||||
function [gamma, maxima, eig_cl] = run_case(Acase, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev)
|
||||
% Solve Lyapunov, compute invariant level, per-halfspace maxima.
|
||||
P = lyap(Acase.', Qbar);
|
||||
Ph = sqrtm(P); Phi = inv(Ph);
|
||||
mu = min(eig(Phi * Qbar * Phi));
|
||||
g_bound = sqrt(B_w.' * P * B_w);
|
||||
c_inv = (2 * w_bar * g_bound / mu)^2;
|
||||
c_entry = delta_entry.' * abs(P) * delta_entry;
|
||||
gamma = max(c_inv, c_entry);
|
||||
Pinv = inv(P);
|
||||
maxima = zeros(size(A_inv, 1), 1);
|
||||
for k = 1:size(A_inv, 1)
|
||||
a = A_inv(k, :).';
|
||||
maxima(k) = sqrt(gamma * (a.' * Pinv * a));
|
||||
end
|
||||
eig_cl = eig(Acase);
|
||||
end
|
||||
|
||||
%% ===== Run both =====
|
||||
[gamma_OL, max_OL, eig_OL] = run_case(A, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev);
|
||||
[gamma_CL, max_CL, eig_CL] = run_case(A_cl, B_w, Qbar, w_bar, delta_entry, A_inv, b_dev);
|
||||
|
||||
%% ===== Report =====
|
||||
fprintf('\n=== Open-loop vs LQR closed-loop Lyapunov barrier ===\n');
|
||||
fprintf(' max Re(eig) OL = %.3e CL = %.3e\n', max(real(eig_OL)), max(real(eig_CL)));
|
||||
fprintf(' min Re(eig) OL = %.3e CL = %.3e\n', min(real(eig_OL)), min(real(eig_CL)));
|
||||
fprintf(' gamma OL = %.3e CL = %.3e (ratio CL/OL = %.3g)\n', ...
|
||||
gamma_OL, gamma_CL, gamma_CL/gamma_OL);
|
||||
|
||||
fprintf('\n Per-halfspace max |a'' dx| on gamma-ellipsoid:\n');
|
||||
fprintf(' %-22s %-10s %-12s %-12s %-12s %-10s\n', ...
|
||||
'halfspace', 'headroom', 'open-loop', 'closed-loop', 'CL - OL', 'ratio');
|
||||
for k = 1:size(A_inv, 1)
|
||||
ratio = max_CL(k) / max_OL(k);
|
||||
fprintf(' %-22s %10.3f %12.3f %12.3f %+12.3f %8.3fx\n', ...
|
||||
comps{k}, b_dev(k), max_OL(k), max_CL(k), max_CL(k) - max_OL(k), ratio);
|
||||
end
|
||||
|
||||
fprintf('\n Interpretation:\n');
|
||||
fprintf(' - CL < OL on a row ==> LQR tightens that halfspace.\n');
|
||||
fprintf(' - CL ~= OL ==> LQR does not help that direction at all.\n');
|
||||
fprintf(' - CL > OL ==> LQR made that direction WORSE for the barrier (!!).\n');
|
||||
@ -1,200 +0,0 @@
|
||||
%% barrier_lyapunov.m — Lyapunov-ellipsoid barrier certificate
|
||||
%
|
||||
% For dx/dt = A_cl x + B_w w with A_cl Hurwitz and ||w||_inf <= w_bar:
|
||||
%
|
||||
% 1. Solve A_cl' P + P A_cl = -Qbar (Qbar > 0, chosen = I).
|
||||
% Then V(x) = x' P x is a Lyapunov function for the undisturbed
|
||||
% system, with dV/dt = -x'x - x'(Qbar-I)x (here Qbar=I gives -x'x).
|
||||
%
|
||||
% 2. Under bounded disturbance:
|
||||
% dV/dt = -x'x + 2 x' P B_w w
|
||||
% <= -||x||^2 + 2 ||P B_w|| w_bar ||x||.
|
||||
% dV/dt <= 0 whenever ||x|| >= 2 ||P B_w|| w_bar.
|
||||
% So the ball B_r := {x : ||x|| <= 2 ||P B_w|| w_bar} contains
|
||||
% the set where V can still grow. Any level set {V <= c} that
|
||||
% contains B_r is forward-invariant.
|
||||
%
|
||||
% 3. Smallest such c: c* = lambda_max(P) * r^2, where r = 2||P B_w||w_bar.
|
||||
%
|
||||
% 4. Safety: the barrier is B(x) = V(x) - gamma, with gamma chosen
|
||||
% large enough to contain X_entry but small enough that the level
|
||||
% set stays inside X_safe. We report whether such a gamma exists.
|
||||
%
|
||||
% This is an ellipsoidal over-approximation, generally much looser than
|
||||
% the box/zonotope reach in reach_operation.m, but it gives a *certificate*
|
||||
% (a closed-form invariant function) rather than just a numerical tube.
|
||||
|
||||
clear; clc;
|
||||
addpath('../plant-model', '../plant-model/controllers');
|
||||
|
||||
plant = pke_params();
|
||||
x_op = pke_initial_conditions(plant);
|
||||
|
||||
%% ===== Build A_cl, B_w =====
|
||||
[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0);
|
||||
|
||||
Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e-2, 1e2, 1]);
|
||||
R_lqr = 1e6;
|
||||
try
|
||||
K = lqr(A, B, Q_lqr, R_lqr);
|
||||
catch
|
||||
[~, ~, K] = icare(A, B, Q_lqr, R_lqr);
|
||||
end
|
||||
A_cl = A - B*K;
|
||||
|
||||
%% ===== Solve Lyapunov equation =====
|
||||
% A_cl' P + P A_cl + Qbar = 0. Qbar shaped to weight T_c heavily so the
|
||||
% resulting ellipsoidal invariant sets are tight in the T_c direction.
|
||||
% Without shaping, isotropic Qbar = I gives ellipsoids stretched along
|
||||
% the slow-precursor directions, making the T_c safety bound useless.
|
||||
Qbar = diag([1, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1e-4, 1, 1e3, 1]);
|
||||
P = lyap(A_cl.', Qbar);
|
||||
assert(all(eig(P) > 0), 'P not positive definite');
|
||||
|
||||
%% ===== Safety spec (used by sweep and final check) =====
|
||||
% Load inv2_holds (conjunction of safety halfspaces) from the predicates
|
||||
% source of truth. Each row k of A_inv defines a halfspace a_k' x <= b_k;
|
||||
% the barrier must bound max(a_k' * dx) over the ellipsoid for each k.
|
||||
addpath('../plant-model');
|
||||
pred = load_predicates(plant);
|
||||
inv2 = pred.mode_invariants.inv2_holds;
|
||||
A_inv = inv2.A_poly;
|
||||
b_inv = inv2.b_poly;
|
||||
comp_names = inv2.components;
|
||||
n_halfspaces = size(A_inv, 1);
|
||||
|
||||
% Convert limits to deviation from x_op:
|
||||
% for halfspace a' x <= b, the deviation-frame bound is a' dx <= b - a' x_op.
|
||||
b_inv_dev = b_inv - A_inv * x_op;
|
||||
|
||||
% Backward-compat scalars for existing prints.
|
||||
e9 = zeros(10, 1); e9(9) = 1;
|
||||
delta_safe_Tc = 5.0;
|
||||
|
||||
%% ===== Disturbance bound =====
|
||||
% |w| <= w_bar where w = Q_sg - Q_nom. Take the same 15% down-load as
|
||||
% reach_operation.m.
|
||||
w_bar = 0.15 * plant.P0;
|
||||
|
||||
% --- Invariant-level computation ---
|
||||
% dV/dt = -x' Qbar x + 2 x' P B_w w.
|
||||
% Taking the worst w = w_bar * sign(x' P B_w), the scalar g = x' P B_w:
|
||||
% dV/dt <= -x' Qbar x + 2 w_bar |g|.
|
||||
% Let u = P^{1/2} x (so V = ||u||^2). Then |g| = |u' P^{-1/2} P B_w|
|
||||
% <= ||u|| * ||P^{-1/2} P B_w|| = sqrt(V) * sqrt(B_w' P B_w).
|
||||
% And x' Qbar x >= lambda_min(P^{-1/2} Qbar P^{-1/2}) * V (call this mu).
|
||||
% So dV/dt <= -mu V + 2 w_bar sqrt(B_w' P B_w) sqrt(V).
|
||||
% dV/dt <= 0 whenever sqrt(V) >= 2 w_bar sqrt(B_w' P B_w) / mu,
|
||||
% i.e. V >= (2 w_bar sqrt(B_w' P B_w) / mu)^2 := c_inv.
|
||||
%
|
||||
% This is much tighter than the isotropic ball bound — it uses the fact
|
||||
% that B_w only pokes one direction of the ellipsoid.
|
||||
P_half = sqrtm(P);
|
||||
P_half_inv = inv(P_half);
|
||||
mu = min(eig(P_half_inv * Qbar * P_half_inv));
|
||||
g_bound = sqrt(B_w.' * P * B_w); % sqrt(B_w' P B_w)
|
||||
c_inv = (2 * w_bar * g_bound / mu)^2;
|
||||
|
||||
fprintf('\n=== Lyapunov barrier certificate ===\n');
|
||||
fprintf(' lambda_min(P) = %.3e\n', min(eig(P)));
|
||||
fprintf(' lambda_max(P) = %.3e\n', max(eig(P)));
|
||||
fprintf(' sqrt(B_w'' P B_w) = %.3e\n', g_bound);
|
||||
fprintf(' mu (Qbar eig on P-metric) = %.3e\n', mu);
|
||||
fprintf(' w_bar (15%% P0) = %.3e W\n', w_bar);
|
||||
fprintf(' c_inv (invariant level) = %.3e\n', c_inv);
|
||||
|
||||
%% ===== Containment of initial set =====
|
||||
% Initial set: box around x_op with halfwidth delta_entry (matches reach_operation).
|
||||
delta_entry = [0.01 * x_op(1);
|
||||
0.001 * abs(x_op(2:7));
|
||||
0.1; 0.1; 0.1];
|
||||
|
||||
% Worst-case V over the initial box: max x'Px = sum over all 2^n corners.
|
||||
% For small n we could enumerate, but the sharper bound is
|
||||
% max V(dx) = delta_entry' * |P| * delta_entry
|
||||
% (elementwise abs of P), which is the L1 energy bound.
|
||||
c_entry = delta_entry.' * abs(P) * delta_entry;
|
||||
|
||||
fprintf('\n c_entry (bound on V over initial box) = %.3e\n', c_entry);
|
||||
|
||||
gamma = max(c_entry, c_inv); % barrier level must contain both
|
||||
fprintf(' gamma (barrier level) = %.3e\n', gamma);
|
||||
if gamma == c_entry
|
||||
fprintf(' (initial set drives gamma — invariant piece already inside entry)\n');
|
||||
else
|
||||
fprintf(' (disturbance inflation drives gamma)\n');
|
||||
end
|
||||
|
||||
%% ===== Sweep Qbar(9,9) to find the tightest safe barrier =====
|
||||
% The isotropic Lyapunov bound is conservative because the "slow decay"
|
||||
% direction dominates mu even when T_c is tightly controlled. Sweep the
|
||||
% T_c weight to find a Qbar that yields a sub-5K barrier if one exists
|
||||
% for this LQR design.
|
||||
fprintf('\n=== Sweeping Qbar(T_c) weight ===\n');
|
||||
weights = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6];
|
||||
best_dTc = inf; best_w = NaN; best_gamma = NaN; best_P = [];
|
||||
for wTc = weights
|
||||
Qbar_s = Qbar; Qbar_s(9,9) = wTc;
|
||||
try
|
||||
Ps = lyap(A_cl.', Qbar_s);
|
||||
catch
|
||||
continue
|
||||
end
|
||||
if any(eig(Ps) <= 0), continue, end
|
||||
Ph = sqrtm(Ps); Phi = inv(Ph);
|
||||
mu_s = min(eig(Phi * Qbar_s * Phi));
|
||||
g_s = sqrt(B_w.' * Ps * B_w);
|
||||
ci_s = (2 * w_bar * g_s / mu_s)^2;
|
||||
ce_s = delta_entry.' * abs(Ps) * delta_entry;
|
||||
g_s_level = max(ci_s, ce_s);
|
||||
Pinv_s = inv(Ps);
|
||||
dTc_s = sqrt(g_s_level * (e9.' * Pinv_s * e9));
|
||||
fprintf(' Qbar(9,9) = %.0e -> gamma = %.3e, max|dT_c| = %7.3f K\n', ...
|
||||
wTc, g_s_level, dTc_s);
|
||||
if dTc_s < best_dTc
|
||||
best_dTc = dTc_s; best_w = wTc; best_gamma = g_s_level; best_P = Ps;
|
||||
end
|
||||
end
|
||||
fprintf(' Best: Qbar(9,9) = %.0e -> max|dT_c| = %.3f K\n', best_w, best_dTc);
|
||||
if best_dTc <= delta_safe_Tc
|
||||
fprintf(' *** TIGHT BARRIER FOUND: V(x) = x.'' P_best x - gamma ***\n');
|
||||
P = best_P; gamma = best_gamma;
|
||||
end
|
||||
|
||||
%% ===== Safety: does the gamma-level set fit inside X_safe? =====
|
||||
% X_safe = { x : |T_c - T_c0| <= 5 K }, i.e. |e_9.' * dx| <= 5.
|
||||
% Max |e_9.' * dx| over {dx : dx' P dx <= gamma} is sqrt(gamma * e_9' P^-1 e_9).
|
||||
Pinv = inv(P);
|
||||
max_dTc_on_ellipsoid = sqrt(gamma * (e9.' * Pinv * e9));
|
||||
|
||||
fprintf('\n=== Safety check on T_c ===\n');
|
||||
fprintf(' Max |dT_c| on gamma-ellipsoid = %.3f K\n', max_dTc_on_ellipsoid);
|
||||
fprintf(' Safe band = +/- %.1f K\n', delta_safe_Tc);
|
||||
if max_dTc_on_ellipsoid <= delta_safe_Tc
|
||||
fprintf(' BARRIER VALID: V(x) = x.''Px - %.3e certifies T_c safety.\n', gamma);
|
||||
else
|
||||
fprintf(' *** BARRIER TOO LOOSE *** - ellipsoid reach into unsafe region.\n');
|
||||
fprintf(' Try a tighter LQR (bigger Q_Tc or smaller R) or tighter X_entry.\n');
|
||||
end
|
||||
|
||||
%% ===== Per-halfspace barrier check against inv2_holds =====
|
||||
% For each safety halfspace a' dx <= b_dev, the max of a' dx over the
|
||||
% gamma-ellipsoid {dx : dx' P dx <= gamma} is sqrt(gamma * a' P^{-1} a).
|
||||
% Compare to b_dev (the headroom from x_op to the safety limit).
|
||||
fprintf('\n=== Lyapunov barrier vs inv2_holds halfspaces ===\n');
|
||||
Pinv = inv(P);
|
||||
for k = 1:n_halfspaces
|
||||
a = A_inv(k, :).';
|
||||
max_adx = sqrt(gamma * (a.' * Pinv * a));
|
||||
margin = b_inv_dev(k) - max_adx;
|
||||
status = 'OK';
|
||||
if margin < 0, status = '*** BARRIER TOO LOOSE ***'; end
|
||||
fprintf(' [%-20s] headroom = %8.3f | max a''dx = %8.3f | margin = %+8.3f %s\n', ...
|
||||
comp_names{k}, b_inv_dev(k), max_adx, margin, status);
|
||||
end
|
||||
|
||||
save(fullfile('.', 'barrier_lyapunov_result.mat'), ...
|
||||
'P', 'gamma', 'c_entry', 'c_inv', 'w_bar', 'K', 'A_cl', 'delta_entry', ...
|
||||
'max_dTc_on_ellipsoid', 'delta_safe_Tc', 'A_inv', 'b_inv_dev', '-v7');
|
||||
|
||||
fprintf('\nSaved barrier to ./barrier_lyapunov_result.mat\n');
|
||||
@ -1,118 +0,0 @@
|
||||
function pred = load_predicates(plant)
|
||||
% LOAD_PREDICATES Read predicates.json and resolve rhs_expr into numbers.
|
||||
%
|
||||
% Returns:
|
||||
% pred.constants - struct with T_c0, T_cold0, T_f0, T_standby
|
||||
% pred.operational_deadbands - struct of predicates (each: A_poly, b_poly, meaning)
|
||||
% pred.safety_limits - struct of halfspace limits (each: A_poly, b_poly, meaning)
|
||||
% pred.mode_invariants - struct mapping mode-invariant name to
|
||||
% its conjoined (A_poly, b_poly, components)
|
||||
%
|
||||
% Each halfspace is of the form { x : A_poly * x <= b_poly }. Conjunctions
|
||||
% (polytopes) stack rows into A_poly / b_poly.
|
||||
%
|
||||
% Usage:
|
||||
% plant = pke_params();
|
||||
% pred = load_predicates(plant);
|
||||
% inv2 = pred.mode_invariants.inv2_holds;
|
||||
% is_in = all(inv2.A_poly * x <= inv2.b_poly);
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
raw = fileread(fullfile(here, 'predicates.json'));
|
||||
J = jsondecode(raw);
|
||||
|
||||
% --- Constants used in rhs_expr evaluations ---
|
||||
T_c0 = plant.T_c0; %#ok<NASGU>
|
||||
T_f0 = plant.T_f0; %#ok<NASGU>
|
||||
T_cold0 = plant.T_cold0; %#ok<NASGU>
|
||||
T_standby_offset_C = J.derived.T_standby_offset_C;
|
||||
T_standby = T_c0 + T_standby_offset_C; %#ok<NASGU>
|
||||
|
||||
pred.constants = struct( ...
|
||||
'T_c0', plant.T_c0, ...
|
||||
'T_f0', plant.T_f0, ...
|
||||
'T_cold0', plant.T_cold0, ...
|
||||
'T_standby', T_standby, ...
|
||||
'T_standby_offset_C', T_standby_offset_C, ...
|
||||
'T_standby_offset_F', J.derived.T_standby_offset_F, ...
|
||||
't_avg_in_range_halfwidth_C', J.derived.t_avg_in_range_halfwidth_C, ...
|
||||
'p_above_crit_threshold_n', J.derived.p_above_crit_threshold_n, ...
|
||||
'T_fuel_limit_C', J.derived.T_fuel_limit_C, ...
|
||||
'T_c_high_trip_C', J.derived.T_c_high_trip_C, ...
|
||||
'n_high_trip', J.derived.n_high_trip);
|
||||
|
||||
% --- operational_deadbands ---
|
||||
pred.operational_deadbands = parse_group(J.operational_deadbands, ...
|
||||
T_c0, T_f0, T_cold0, T_standby);
|
||||
|
||||
% --- safety_limits ---
|
||||
pred.safety_limits = parse_group(J.safety_limits, ...
|
||||
T_c0, T_f0, T_cold0, T_standby);
|
||||
|
||||
% --- mode_invariants: conjunctions of safety_limits entries ---
|
||||
inv_names = fieldnames(J.mode_invariants);
|
||||
for k = 1:numel(inv_names)
|
||||
name = inv_names{k};
|
||||
if startsWith(name, '_') || startsWith(name, 'x_'), continue, end
|
||||
entry = J.mode_invariants.(name);
|
||||
if ~isstruct(entry) || ~isfield(entry, 'conjunction_of'), continue, end
|
||||
components = entry.conjunction_of;
|
||||
if ischar(components), components = {components}; end
|
||||
A_all = [];
|
||||
b_all = [];
|
||||
for i = 1:numel(components)
|
||||
comp = components{i};
|
||||
A_all = [A_all; pred.safety_limits.(comp).A_poly]; %#ok<AGROW>
|
||||
b_all = [b_all; pred.safety_limits.(comp).b_poly]; %#ok<AGROW>
|
||||
end
|
||||
pred.mode_invariants.(name).A_poly = A_all;
|
||||
pred.mode_invariants.(name).b_poly = b_all;
|
||||
pred.mode_invariants.(name).meaning = entry.meaning;
|
||||
pred.mode_invariants.(name).components = components;
|
||||
end
|
||||
end
|
||||
|
||||
function group_out = parse_group(group_in, T_c0, T_f0, T_cold0, T_standby)
|
||||
names = fieldnames(group_in);
|
||||
group_out = struct();
|
||||
for k = 1:numel(names)
|
||||
name = names{k};
|
||||
% MATLAB jsondecode renames "_comment" -> "x_comment" and similar
|
||||
if startsWith(name, '_') || startsWith(name, 'x_'), continue, end
|
||||
entry = group_in.(name);
|
||||
if ~isstruct(entry) || ~isfield(entry, 'halfspaces'), continue, end
|
||||
hs_list = entry.halfspaces;
|
||||
if iscell(hs_list)
|
||||
n_hs = numel(hs_list);
|
||||
get_hs = @(i) hs_list{i};
|
||||
else
|
||||
n_hs = numel(hs_list);
|
||||
get_hs = @(i) hs_list(i);
|
||||
end
|
||||
A_poly = zeros(n_hs, 10);
|
||||
b_poly = zeros(n_hs, 1);
|
||||
for i = 1:n_hs
|
||||
hs = get_hs(i);
|
||||
if isfield(hs, 'row')
|
||||
% Multi-coefficient halfspace: hs.row is a k x 2 matrix
|
||||
% where each row is [state_index, coeff] (jsondecode from
|
||||
% JSON array-of-arrays).
|
||||
coeffs = hs.row;
|
||||
for r = 1:size(coeffs, 1)
|
||||
A_poly(i, coeffs(r, 1)) = coeffs(r, 2);
|
||||
end
|
||||
else
|
||||
% Single-coefficient halfspace: hs.state_index + hs.coeff.
|
||||
A_poly(i, hs.state_index) = hs.coeff;
|
||||
end
|
||||
b_poly(i) = evalin_context(hs.rhs_expr, T_c0, T_f0, T_cold0, T_standby);
|
||||
end
|
||||
group_out.(name).A_poly = A_poly;
|
||||
group_out.(name).b_poly = b_poly;
|
||||
group_out.(name).meaning = entry.meaning;
|
||||
end
|
||||
end
|
||||
|
||||
function val = evalin_context(expr, T_c0, T_f0, T_cold0, T_standby) %#ok<INUSD>
|
||||
val = eval(expr);
|
||||
end
|
||||
@ -1,93 +0,0 @@
|
||||
function [T, R_lo, R_hi, center] = reach_linear(A_cl, B_w, x0_center, x0_halfwidth, w_lo, w_hi, tspan, dt)
|
||||
% REACH_LINEAR Interval reach set for dx/dt = A_cl*x + B_w*w, w in [w_lo, w_hi].
|
||||
%
|
||||
% Hand-rolled zonotope propagator specialized to the case where:
|
||||
% - The initial set is an axis-aligned box around x0_center with
|
||||
% halfwidths x0_halfwidth.
|
||||
% - The disturbance w is a scalar interval.
|
||||
%
|
||||
% The reach set at time t is:
|
||||
% X(t) = {expm(A_cl*t) * x0 + integral_0^t expm(A_cl*(t-s))*B_w*w(s) ds :
|
||||
% x0 in X0, w(s) in W}.
|
||||
%
|
||||
% Decompose:
|
||||
% X(t) = expm(A_cl*t) * X0 (+) S_w(t),
|
||||
% where S_w(t) is the reachable contribution of the disturbance.
|
||||
% Both parts are zonotopes; we over-approximate each by a box at every
|
||||
% time step via the matrix/integral absolute sum (interval hull).
|
||||
%
|
||||
% Inputs:
|
||||
% A_cl - n x n closed-loop state matrix (Hurwitz assumed)
|
||||
% B_w - n x 1 disturbance vector
|
||||
% x0_center - n x 1 center of initial set
|
||||
% x0_halfwidth - n x 1 halfwidths of initial set (>=0)
|
||||
% w_lo, w_hi - scalar disturbance bounds
|
||||
% tspan - [t0, tf]
|
||||
% dt - time step for exporting the reach tube (s)
|
||||
%
|
||||
% Outputs:
|
||||
% T - M x 1 time grid
|
||||
% R_lo - n x M lower bounds of reach set
|
||||
% R_hi - n x M upper bounds of reach set
|
||||
% center- n x M center trajectory (nominal, w = w_mid)
|
||||
|
||||
n = size(A_cl, 1);
|
||||
assert(size(A_cl,2) == n, 'A_cl must be square');
|
||||
assert(numel(x0_center) == n && numel(x0_halfwidth) == n, 'x0 dim mismatch');
|
||||
x0_center = x0_center(:);
|
||||
x0_halfwidth = x0_halfwidth(:);
|
||||
B_w = B_w(:);
|
||||
|
||||
w_mid = 0.5 * (w_lo + w_hi);
|
||||
w_half = 0.5 * (w_hi - w_lo);
|
||||
|
||||
T = (tspan(1):dt:tspan(2)).';
|
||||
M = numel(T);
|
||||
|
||||
R_lo = zeros(n, M);
|
||||
R_hi = zeros(n, M);
|
||||
center = zeros(n, M);
|
||||
|
||||
% Disturbance contribution:
|
||||
% S_center(t) = int_0^t expm(A*(t-s))*B_w*w_mid ds (signed, centerline)
|
||||
% S_half(t) = int_0^t |expm(A*(t-s))*B_w| * w_half ds (elementwise, halfwidth)
|
||||
%
|
||||
% The halfwidth uses elementwise |.| because under an axis-aligned
|
||||
% box over-approximation, |linear-map| contributes absolute values.
|
||||
% This is the interval-arithmetic analog of propagating a zonotope.
|
||||
S_center = zeros(n, 1);
|
||||
S_half = zeros(n, 1);
|
||||
|
||||
x_center_now = x0_center;
|
||||
A_step = expm(A_cl * dt); % state transition per step
|
||||
A_abs_step = abs(A_step); % for interval-halfwidth
|
||||
|
||||
% Integrated input gain per step:
|
||||
% G_step = int_0^dt expm(A*s)*B_w ds
|
||||
% Augmented-matrix trick: upper-right block of expm([A B_w; 0 0]*dt).
|
||||
M_aug = expm([A_cl, B_w; zeros(1, n+1)] * dt);
|
||||
G_step = M_aug(1:n, n+1); % n x 1
|
||||
G_abs_step = abs(G_step);
|
||||
|
||||
halfwidth_now = x0_halfwidth;
|
||||
|
||||
% Initial state
|
||||
center(:, 1) = x_center_now + S_center;
|
||||
R_lo(:, 1) = center(:, 1) - halfwidth_now - S_half;
|
||||
R_hi(:, 1) = center(:, 1) + halfwidth_now + S_half;
|
||||
|
||||
for k = 2:M
|
||||
% State piece
|
||||
x_center_now = A_step * x_center_now;
|
||||
halfwidth_now = A_abs_step * halfwidth_now;
|
||||
|
||||
% Disturbance piece
|
||||
S_center = A_step * S_center + G_step * w_mid;
|
||||
S_half = A_abs_step * S_half + G_abs_step * w_half;
|
||||
|
||||
center(:, k) = x_center_now + S_center;
|
||||
R_lo(:, k) = center(:, k) - halfwidth_now - S_half;
|
||||
R_hi(:, k) = center(:, k) + halfwidth_now + S_half;
|
||||
end
|
||||
|
||||
end
|
||||
@ -1,203 +0,0 @@
|
||||
%% reach_operation.m — linear reach set for operation mode (LQR closed-loop)
|
||||
%
|
||||
% *** SOUNDNESS STATUS: APPROXIMATE, NOT SOUND. ***
|
||||
%
|
||||
% This file computes a reach tube for the *linearized* closed-loop
|
||||
% system (A_cl = A - BK around x_op) under bounded Q_sg. The tube
|
||||
% itself is a sound over-approximation of the LINEAR model's reach
|
||||
% set — it uses conservative box hulls and elementwise-absolute-value
|
||||
% matrix propagation. But the LINEAR MODEL is only an approximation
|
||||
% of the real nonlinear plant (pke_th_rhs.m), so the result is not a
|
||||
% sound reach tube for the actual plant.
|
||||
%
|
||||
% To upgrade to a sound result, pick one:
|
||||
% (a) Nonlinear reach directly (CORA's nonlinearSys, JuliaReach's
|
||||
% BlackBoxContinuousSystem). Expensive, honest.
|
||||
% (b) Linear reach + Taylor-remainder inflation: compute an upper
|
||||
% bound on ||f_nonlinear(x,u) - (A x + B u)|| over the reach
|
||||
% set and inflate the linear tube by that bound. Cheaper,
|
||||
% requires a Hessian-norm estimate.
|
||||
% Tracked as a thesis-blocking todo. For now, the 5-orders-of-margin
|
||||
% buffer (|dT_c| ~ 0.03 K vs safety band 5 K) gives us a lot of room
|
||||
% to absorb linearization error, but that's not a proof.
|
||||
%
|
||||
% Compute a reach-tube over-approximation starting from a box around
|
||||
% x_op, under LQR feedback, with Q_sg in a specified interval. Check
|
||||
% that T_avg stays inside the t_avg_in_range predicate for all t in
|
||||
% [0, T_final].
|
||||
%
|
||||
% This is the *continuous-mode obligation* for q_operation:
|
||||
% X_entry := { x : |x - x_op| <= delta_entry }
|
||||
% W := [Q_min, Q_max]
|
||||
% X_safe := { x : |T_c - T_c0| <= delta_safe }
|
||||
% Obligation: ReachTube(X_entry, W, [0, T_final]) subset X_safe.
|
||||
%
|
||||
% If this passes, we've discharged the Thrust-3 verification for one
|
||||
% continuous mode at the level the thesis calls for.
|
||||
|
||||
clear; clc; close all;
|
||||
addpath('../plant-model', '../plant-model/controllers');
|
||||
|
||||
plant = pke_params();
|
||||
x_op = pke_initial_conditions(plant);
|
||||
pred = load_predicates(plant); % single source of truth for predicate bands
|
||||
|
||||
%% ===== Closed-loop linearization =====
|
||||
[A, B, B_w, ~, ~, ~] = pke_linearize(plant, x_op, 0, plant.P0);
|
||||
|
||||
% LQR gain using the same Q, R as ctrl_operation_lqr.m. Keep these in
|
||||
% sync if you retune there.
|
||||
Q_lqr = diag([1, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-2, 1e2, 1]);
|
||||
R_lqr = 1e6;
|
||||
try
|
||||
K = lqr(A, B, Q_lqr, R_lqr);
|
||||
catch
|
||||
[~, ~, K] = icare(A, B, Q_lqr, R_lqr);
|
||||
end
|
||||
|
||||
A_cl = A - B*K;
|
||||
|
||||
fprintf('\n=== Closed-loop spectrum (A - BK) ===\n');
|
||||
eigs_cl = eig(A_cl);
|
||||
fprintf(' max Re(eig) = %.3e\n', max(real(eigs_cl)));
|
||||
fprintf(' min Re(eig) = %.3e\n', min(real(eigs_cl)));
|
||||
assert(all(real(eigs_cl) < -1e-8), 'A_cl not Hurwitz — gain tuning issue');
|
||||
|
||||
%% ===== Define sets =====
|
||||
% X_entry: 1% box on n, 0.1% boxes on precursors (their magnitudes are
|
||||
% huge due to 1/Lambda), 0.1 K on fuel, 0.1 K on coolant, 0.1 K on cold leg.
|
||||
% This represents "we entered operation mode near the steady state".
|
||||
delta_entry = [0.01 * x_op(1); % n
|
||||
0.001 * abs(x_op(2:7)); % C1..C6
|
||||
0.1; % T_f [C]
|
||||
0.1; % T_c [C]
|
||||
0.1]; % T_cold [C]
|
||||
|
||||
% Disturbance: Q_sg = [0.85*P0, 1.00*P0] -> this captures up to a 15%
|
||||
% down-step load demand (realistic load-follow envelope).
|
||||
Q_nom = plant.P0;
|
||||
Q_min = 0.85 * plant.P0;
|
||||
Q_max = 1.00 * plant.P0;
|
||||
dQ_lo = Q_min - Q_nom; % -0.15 * P0
|
||||
dQ_hi = Q_max - Q_nom; % 0
|
||||
|
||||
% X_safe is inv2_holds — the operation-mode safety envelope. Each row of
|
||||
% inv2.A_poly is a hard safety limit (fuel centerline, T_c high/low trip,
|
||||
% n high/low, cold-leg subcooling). We check reach-tube containment
|
||||
% halfspace-by-halfspace so failure modes are attributable.
|
||||
inv2 = pred.mode_invariants.inv2_holds;
|
||||
% Keep the old +-halfwidth report too, for continuity with the
|
||||
% operational-deadband framing.
|
||||
delta_safe_Tc = pred.constants.t_avg_in_range_halfwidth_C; % [C]
|
||||
|
||||
%% ===== Reach set =====
|
||||
% Propagate in deviation coordinates: dx = x - x_op.
|
||||
tspan = [0, 600];
|
||||
dt = 0.5;
|
||||
|
||||
[T, R_lo, R_hi, C] = reach_linear(A_cl, B_w, zeros(10,1), delta_entry, dQ_lo, dQ_hi, tspan, dt);
|
||||
|
||||
% Translate back to absolute coordinates for reporting
|
||||
Xabs_lo = R_lo + x_op;
|
||||
Xabs_hi = R_hi + x_op;
|
||||
Cabs = C + x_op;
|
||||
|
||||
%% ===== Safety check =====
|
||||
T_c_lo = Xabs_lo(9, :);
|
||||
T_c_hi = Xabs_hi(9, :);
|
||||
|
||||
violation_mask = (T_c_hi > plant.T_c0 + delta_safe_Tc) | ...
|
||||
(T_c_lo < plant.T_c0 - delta_safe_Tc);
|
||||
fprintf('\n=== Operation-mode reach-set safety ===\n');
|
||||
fprintf(' Horizon = [%g, %g] s\n', tspan(1), tspan(2));
|
||||
fprintf(' Entry box T_c [C] = [%.3f, %.3f] (x_op +/- %.1f C)\n', ...
|
||||
x_op(9) - delta_entry(9), x_op(9) + delta_entry(9), delta_entry(9));
|
||||
fprintf(' Disturbance Q_sg = [%.3f, %.3f] MW\n', Q_min/1e6, Q_max/1e6);
|
||||
fprintf(' Safe band on T_c = x_op(T_c0) +/- %.1f C -> [%.3f, %.3f]\n', ...
|
||||
delta_safe_Tc, plant.T_c0 - delta_safe_Tc, plant.T_c0 + delta_safe_Tc);
|
||||
fprintf(' Reach T_c envelope = [%.3f, %.3f]\n', min(T_c_lo), max(T_c_hi));
|
||||
if any(violation_mask)
|
||||
t_first = T(find(violation_mask, 1));
|
||||
fprintf(' *** SAFETY VIOLATED at t = %.2f s ***\n', t_first);
|
||||
else
|
||||
fprintf(' OK: reach set stays inside the safe band.\n');
|
||||
end
|
||||
|
||||
%% Hard safety-limit check (inv2_holds halfspace-by-halfspace)
|
||||
% For each row a_k of inv2.A_poly with threshold b_k, check whether
|
||||
% max over reach tube of a_k * x stays <= b_k. The reach tube upper
|
||||
% envelope is Xabs_hi; lower envelope is Xabs_lo. We evaluate
|
||||
% max(a_k * x) using Xabs_hi where a_k > 0, Xabs_lo where a_k < 0.
|
||||
fprintf('\n=== Operation-mode reach vs inv2_holds safety limits ===\n');
|
||||
A_inv = inv2.A_poly; b_inv = inv2.b_poly;
|
||||
comps = inv2.components;
|
||||
for k = 1:size(A_inv, 1)
|
||||
a = A_inv(k, :).';
|
||||
% envelope maximum of a' * x across the reach tube
|
||||
x_envelope = Xabs_hi .* (a > 0) + Xabs_lo .* (a < 0); % 10 x M
|
||||
max_ax = max(a.' * x_envelope);
|
||||
margin = b_inv(k) - max_ax;
|
||||
status = 'OK';
|
||||
if margin < 0, status = '*** VIOLATED ***'; end
|
||||
fprintf(' [%s] a''x <= %.3f | max a''x = %.3f | margin = %+.3f %s\n', ...
|
||||
comps{k}, b_inv(k), max_ax, margin, status);
|
||||
end
|
||||
|
||||
%% Per-state reach-set growth diagnostic (final time vs initial)
|
||||
state_names = {'n','C1','C2','C3','C4','C5','C6','T_f','T_c','T_cold'};
|
||||
fprintf('\n=== Reach-set width at t=0 vs t=T_final ===\n');
|
||||
fprintf(' %-7s %-14s %-14s %-8s\n', 'state', 'init halfwidth', 'final halfwidth', 'ratio');
|
||||
for i = 1:10
|
||||
hi = 0.5 * (R_hi(i, 1) - R_lo(i, 1));
|
||||
hf = 0.5 * (R_hi(i, end) - R_lo(i, end));
|
||||
fprintf(' %-7s %-14.4e %-14.4e %-8.2f\n', state_names{i}, hi, hf, hf/max(hi,eps));
|
||||
end
|
||||
|
||||
%% ===== Plots =====
|
||||
figdir = fullfile('..', 'docs', 'figures');
|
||||
if ~exist(figdir, 'dir'), mkdir(figdir); end
|
||||
CtoF = @(T) T*9/5 + 32;
|
||||
|
||||
% Two-panel plot: wide view with safety band, zoom view showing actual tube.
|
||||
figure('Position', [100 80 1400 500], 'Name', 'Reach tube: T_c');
|
||||
|
||||
subplot(1,2,1);
|
||||
fill([T; flipud(T)], CtoF([T_c_hi.'; flipud(T_c_lo.')]), [1.0 0.85 0.85], ...
|
||||
'EdgeColor', 'none'); hold on;
|
||||
plot(T, CtoF(Cabs(9, :)), 'r-', 'LineWidth', 1.2);
|
||||
yline(CtoF(plant.T_c0 + delta_safe_Tc), 'k--', 'LineWidth', 1.0);
|
||||
yline(CtoF(plant.T_c0 - delta_safe_Tc), 'k--', 'LineWidth', 1.0);
|
||||
yline(CtoF(plant.T_c0), 'k:', 'LineWidth', 1.0);
|
||||
xlabel('Time [s]'); ylabel('T_{avg} [F]'); grid on;
|
||||
title('Safety-band view');
|
||||
legend('reach tube', 'nominal', 'safety +/- 5 C', 'Location', 'best');
|
||||
|
||||
subplot(1,2,2);
|
||||
Tc_dev_lo = T_c_lo.' - plant.T_c0; % M x 1, deviation in K
|
||||
Tc_dev_hi = T_c_hi.' - plant.T_c0;
|
||||
fill([T; flipud(T)], [Tc_dev_hi; flipud(Tc_dev_lo)], [1.0 0.85 0.85], ...
|
||||
'EdgeColor', 'none'); hold on;
|
||||
plot(T, Cabs(9, :).' - plant.T_c0, 'r-', 'LineWidth', 1.2);
|
||||
yline(0, 'k:', 'LineWidth', 1.0);
|
||||
xlabel('Time [s]'); ylabel('T_{avg} - T_{c0} [K]'); grid on;
|
||||
max_dev = max(abs([Tc_dev_lo; Tc_dev_hi]));
|
||||
title(sprintf('Zoomed: max |dT_c| = %.3e K', max_dev));
|
||||
|
||||
sgtitle(sprintf('Operation-mode reach tube, LQR, Q_{sg} in [%.0f%%, %.0f%%] P_0', ...
|
||||
100*Q_min/Q_nom, 100*Q_max/Q_nom));
|
||||
exportgraphics(gcf, fullfile(figdir, 'reach_operation_Tc.png'), 'Resolution', 150);
|
||||
|
||||
figure('Position', [100 80 1100 500], 'Name', 'Reach tube: n');
|
||||
fill([T; flipud(T)], [R_hi(1,:).' + x_op(1); flipud(R_lo(1,:).' + x_op(1))], ...
|
||||
[0.85 0.85 1.0], 'EdgeColor', 'none'); hold on;
|
||||
plot(T, Cabs(1, :), 'b-', 'LineWidth', 1.2);
|
||||
xlabel('Time [s]'); ylabel('n'); grid on;
|
||||
title('Operation mode reach tube on normalized power');
|
||||
legend('reach tube', 'nominal', 'Location', 'best');
|
||||
exportgraphics(gcf, fullfile(figdir, 'reach_operation_n.png'), 'Resolution', 150);
|
||||
|
||||
save(fullfile('.', 'reach_operation_result.mat'), ...
|
||||
'T', 'R_lo', 'R_hi', 'C', 'Xabs_lo', 'Xabs_hi', 'Cabs', ...
|
||||
'K', 'A_cl', 'x_op', 'delta_entry', 'Q_min', 'Q_max', 'delta_safe_Tc', '-v7');
|
||||
|
||||
fprintf('\nSaved reach result to ./reach_operation_result.mat\n');
|
||||
Loading…
x
Reference in New Issue
Block a user