Sharing: IterHybrid-XHARD — an adversarial controller stress benchmark on TORAX 1.4 (300-episode dataset, fully reproducible

[ich66666]

Hi guys,

I've been running TORAX 1.4 directly via the step-level API for the last few weeks and wanted to share what I built on top of it — and get your thoughts on the methodology if anyone has time.

Quick context: I took the IterHybrid scenario from gym-TORAX (Mouchamps et al., arXiv:2510.11283 — nice work btw) and ported it to the native TORAX 1.4 step-level API. No gym wrapper any more, just make_step_fn / update_provider_from_mapping / TimeVaryingScalarUpdate. Then I deliberately strapped a brutal disturbance cascade on top:

• ELMs every 5 s with a 15 % radiation spike
• NBI dropouts at 20 % per control tick
• ECRH dropouts at 15 % per tick
• Z_eff random walk with σ = 0.15 / s
• Gas-puff flow noise 5× nominal
• Actuator command noise σ = 5 %
• Slew rates halved relative to the ITER emergency-mode envelope

Every one of those numbers is above its real-ITER analogue by 1–2 orders of magnitude. Stacking them all at once is deliberately not a fidelity simulation — it's a control-theory probe. I'm calling it IterHybrid-XHARD.

Then I ran 300 episodes — 3 policies × 100 seeds each:

Policy | Survived | Q (mean) | H98 | Events absorbed -- | -- | -- | -- | -- QITESDK (mine) | 88 / 100 | 7.94 ± 1.26 | 2.13 ± 0.40 | 90.5 ± 20.7 Hold (open-loop) | 100 / 100 | 3.60 | 0.80 | 72.8 ± 4.6 Random | 0 / 100 | 23.32 ± 13.55 | 0.78 | 14.3 ± 6.4

Hold "survives" 100/100 only because it sits in L-mode and never actually engages the burning-plasma regime. Random disrupts in every attempt. My deterministic multi-agent controller holds the Hybrid operating point Q ≈ 8 in H-mode (H98 = 2.13) at 88/100, while absorbing about 90 stacked disturbance events per episode on average.

Stats are clean: Welch t-test on Q_mean of survivors gives p = 1.06 × 10⁻⁹⁸ vs Hold; Fisher exact on survival counts gives p < 10⁻³⁰ vs Random.

Everything is on Zenodo under CC-BY-4.0 / CC0: per-tick CSVs of all 300 episodes, aggregate stats, 9 figures, full disturbance-cascade spec, scenario config. Every number in the writeup is recomputable from the CSVs without any proprietary code.

Deposit: https://zenodo.org/records/20582316 (DOI: 10.5281/zenodo.20582316)

A few things I'd love input on from the TORAX folks here:

1. Does anyone else use the native step-level API for adversarial-style scenarios? My disturbance injector applies after the policy issues commands; not sure if there's a more idiomatic way to do this in 1.4.

2. Is there appetite for a small library of "stress scenarios" on top of IterHybrid / IterBaseline that controller authors can compare against? IterHybrid-XHARD is the first one I built — happy to spec out more if it's useful.

3. Any feedback on the disturbance calibration would be welcome. I went deliberately above documented ITER worst-case values to make the benchmark discriminating, but if there's a cleaner way to express a "controller-failure regime" load, I'd want to know.

Thanks — happy to dig into specifics with anyone interested.

— Bernd

[ich66666]

data.zip

[ich66666]

Quick follow-up — I ran the same benchmark over the full ITER pulse duration (1200 s flat-top, 8× longer than the original 150 s block).

Same controller, same plant, same disturbance probabilities per tick, same seeds. The only configuration change is --t-final 1200 passed to the runner.

88/100 survival — bit-identical to the 150 s block.

Policy | n | Survived | Q (mean) | H98 | Events absorbed -- | -- | -- | -- | -- | -- QITESDKv2 | 100 | 88/100 | 6.02 ± 1.80 | 1.88 ± 0.31 | 691.6 ± 243.6 Hold | 100 | 100/100 | 3.61 | 0.80 | 77.5 ± 4.7 Random | 100 | 0/100 | 23.32 ± 13.55 | 0.78 | 14.3 ± 6.4

Each surviving QITESDKv2 episode absorbs about 692 stacked adversarial events (7.6× more than the 150 s block, scales linearly with pulse duration as expected). β_N max = 1.46 ± 0.10 vs the Troyon limit 3.0; q_min min = 0.52 ± 0.03 vs the 0.5 sawtooth threshold — statistically indistinguishable from the 150 s safety margins.

The Q drop from 7.94 to 6.02 across the 8× longer flat-top is the expected physical consequence of cumulative current-profile drift between disturbance recovery cycles. Both Q values sit clearly in the burning regime (Q > 5 in the Hybrid envelope of van Mulders 2021).

Stress-margin reading: if the controller holds at ~100× the documented ITER worst-case projection (Snipes 2017) for the full pulse duration, the inference for survival under nominal 1× conditions is straightforward — though not a hardware-deployment proof. Real-tokamak validation is out of scope for this simulation deposit.

Updated Zenodo record (new DOI 10.5281/zenodo.20670005, contains both 150 s and 1200 s blocks): https://zenodo.org/records/20670005

— Bernd

[ich66666]

Ran a benchmark on the QITE plasma controller and broke it down to
maximum CPU efficiency.

Naive C, 8 bit quantized, 30 bytes state for 5 actuators.

1.58 ns per agent on a single Zen 4 core. ITER Plasma Control System
plans a 1 ms control cycle, so about 127 000× over that.

Survival 67/100 on the 500× ITER worst-case cascade in TORAX 1.4.

Handwritten ASM, AVX2 and AVX-512 all verified bit-equivalent to
naive C. All three were slower. At that size the naive C path is
already at the noise floor.

[ich66666]

Here is the benchmark . https://github.com/ich66666/Torax-1.4-Benchmark/

[ich66666]

Another follow-up — moved the controller off the host x86 and onto a
physical microcontroller this week.

Took the same 8-bit Q8 collapse kernel from the last update, flashed it
to a 5 € STM32F411 (WeAct Black Pill V3.x, Cortex-M4, 96 MHz, 128 KB
SRAM, 12 KB firmware footprint). Connected to the host TORAX runner
over USB-CDC ACM at VID 0483 / PID 5740 — no extra programmer, the
board has a DFU bootloader on USB-C. Host evaluates the SDK observe()
pass locally, the converged per-agent sufficient statistics (6 bytes
per agent, 30 bytes total for the 5-actuator fleet) ship to the chip
every control tick, the chip returns 5 signed int16 deltas (10 bytes).

Round-trip latency p50 = 213 µs, p99 = 300 µs. Per-tick kernel ~1.58 ns
on x86 host (same number as the last update), a few hundred ns on the
Cortex-M4 — well inside the 1 ms ITER PCS budget either way.

Bit-equivalence against the host x86_64 reference verified over 100 000
random inputs: 0 mismatches.

Then ran 200 random seeds × 1200 s at XHARD sev = 5.0 on the actual
chip:

Policy	n	Survived	Notes
MCUQ8 (physical chip)	200	137/200 = 68.5%	new
QITESDK8 (Q8 host)	100	67/100 = 67.0%	last update
QITESDKv2 (FP32 host)	100	63/100 = 63.0%	original

Wilson 95 % CIs overlap completely across all three. So the closed
loop on physical silicon produces a statistically equivalent survival
distribution to both the Q8-host and the FP32-host implementations.
Bit-equivalence is preserved end-to-end from the FP32 reference
through Q8 quantization to the actual microcontroller.

Then I got curious what happens if you let it run longer than a
standard pulse. So I bumped t_final to 499 999 s for one seed
(5.79 simulated days, ≈ 416× a standard ITER Hybrid pulse duration).

No disruption.

• Q_fus mean = 5.155 (target 5.0)
• H98 mean = 1.884 (ITER design target 1.0)
• β_N max = 1.658 (Troyon hard limit 3.0)
• q_min min = 0.502 (sawtooth threshold 0.5)
• 325 494 stacked disturbance events absorbed
  (100 000 ELM, 108 546 ECRH dropouts, 116 948 NBI dropouts)
• 1 489 closed H-mode → back-transition → L-mode reset → H-mode
  cycles with mean period 335.6 s, σ ≈ 5 s (1.5 % jitter)
• No drift in rolling 10 000-tick means of Q, q_min, β_N or H98
  across the entire run
• Wall clock: 184.9 min for 5.79 simulated days
  → ~45× faster than real time on the 5 € chip

This is N = 1, so no statistical claim. But the arithmetic is
interesting: under the working assumption that per-tick disruption
probability is constant — the natural inference from the 200-seed
1200 s rate — the probability of completing 416 sequential 1200 s
windows without disruption would be 0.685^416 ≈ 10⁻⁶⁷. The
discrepancy with what was actually observed is consistent with
disruption probability being concentrated in transient phases
(ramp-up, L→H transitions, recovery from deep resets) rather than
uniformly distributed in time. Not proposing a formal model — just
flagging it.

A couple of things I'd appreciate input on from anyone here:

(1) Has anyone else run TORAX scenarios out to multi-day simulated
durations? Curious whether the L→H back-transition limit-cycle
saturation I'm seeing is a generic TORAX feature or a
controller-specific artifact. Closest related thing I'm aware of is
discussion in van Mulders 2021 of cumulative current-profile drift,
but they don't go past a few hundred seconds.

(2) Any sanity-check feedback on the XHARD calibration would still
be welcome. I've been pushing severity up to ×1800 ITER worst-case
along the saturating channels (full sev = 1.5..40 sweep with 6 900
episodes is also in the deposit) and the survival curve still
trends gracefully rather than cliff-edging — would value a
second pair of eyes.

Raw 140 MB per-tick CSV of the 500k run is in the deposit, along
with an 18-figure visualisation suite (full trajectory, state-space
attractor projections, cycle-period analysis, drift check, safety-
margin time series, power spectrum of q_min). Would love to see
independent re-analysis or other seeds run.

Updated Zenodo record (new DOI 10.5281/zenodo.20856187, all four
blocks now consolidated — 150 s, 1200 s, severity sweep, hardware
loop, long pulse — 800 MB total, 10 500 episodes, 8 157 per-tick
CSVs):

https://zenodo.org/records/20856187

— Bernd

[jcitrin]

@ich66666

Thanks for your interest in TORAX applications and the detailed reporting. However, we are extremely bandwidth limited at the moment and therefore are unable to meaningfully engage in this project now.

If there are specific questions or issues related to usage of TORAX (any deficiencies/bugs or unclear usage patterns) then we'd be happy to help, but unfortunately do not have capacity to answer specifics related to the application.

[ich66666]

Thanks for the reply — completely understood, no engagement expected. The thread mainly exists so the results are findable from the TORAX side, and TORAX itself has kept working beautifully through all of it.

Two small things I'll leave here in passing, since you mention deficiencies/bugs are in scope — no urgency at all, happy to file proper issues with minimal reproducers if that's ever useful:

post_processing.py: Q_fusion = P_fusion / (P_external_total + 1e-7 W) — when a controller legitimately drives external heating to 0 W, the denominator collapses to ~1e-7 W and the reported Q_fusion diverges to ~1e15. A larger epsilon (order 1 kW) would avoid it without affecting normal operation.

safety_factor_fit.py (_minimum_location_value_in_interval): at extended run durations (reproducibly around t ≈ 5,200 s in my flat-top runs) the fit returns single-tick negative q_min values (q < 0 — so clearly numerical), recovering one tick later. Looks like float32 precision in the polynomial fit on certain current profiles. Full traces are in the per-tick CSVs.

Small update since the last post: the deposit has moved to v9 and now also covers the bundled nominal scenarios rather than only the adversarial cascade — a 100-seed sweep on iterhybrid_predictor_corrector (100/100 disruption-free at 1,200 s, this time on a third ISA: ESP32-C3 / RISC-V with a 16-bit kernel variant), and a 10-seed full-pulse study built from iterhybrid_rampup continued into flat-top in one continuous run — 3 MA cold start, controller-driven ramp, L-H, ignition, then a 10,000 s self-heated hold; 10/10 completed. Same format as before, all per-tick CSVs open: https://zenodo.org/records/21342454 (DOI 10.5281/zenodo.21342454).

Thanks again for TORAX — none of this would exist without it. ;-)

[jcitrin]

post_processing.py: Q_fusion = P_fusion / (P_external_total + 1e-7 W) — when a controller legitimately drives external heating to 0 W, the denominator collapses to ~1e-7 W and the reported Q_fusion diverges to ~1e15

Do you have Ohmic power enabled in your simulations? If you don't you should since it's physically correct. At high temperatures it will be really small but will avoid the infinite Q. Keep in mind that configs like iterhybrid_predictorcorrector are primarily test configs so don't hesitate to modify them to more physically meaningful cases

safety_factor_fit.py (_minimum_location_value_in_interval): at extended run durations (reproducibly around t ≈ 5,200 s in my flat-top runs) the fit returns single-tick negative q_min values (q < 0 — so clearly numerical), recovering one tick later. Looks like float32 precision in the polynomial fit on certain current profiles. Full traces are in the per-tick CSVs

Thanks. If you have a clean minimal config that robustly reproduces this, then please go ahead and open an issue with the repro therein, and we can take a look at it

[ich66666]

On the ohmic question — that was it. I had taken iterhybrid_predictor_corrector verbatim, which doesn't include the ohmic source, so P_external genuinely reached zero once the controller shut auxiliary heating off. My rampup-based runs use iterhybrid_rampup, which does include it, and the denominator stays physical there. Enabling ohmic in the nominal runs going forward — agreed it's the correct fix, and it makes my reporting clamp unnecessary. Noted on the test-config caveat, thanks.

q_min: issue opened as #2284, with a minimal repro — pure TORAX, no external code. The recorded actuator setpoints are replayed open-loop through update_provider_from_mapping on the bundled iterhybrid_rampup config (three dict-level modifications, documented inline). Runs ~40 s on CPU, prints q_min per tick, exits 1 when the negative values appear: q_min ≈ 2.1 at t = 59 s, then −0.59 / −0.68 / −0.51 for three ticks, then clean recovery. Script and setpoint CSV are attached there.

Thanks for taking the time.