On the Computational Limits of AI4S-RL : A Unified $\varepsilon$-$N$ Analysis
Authors: Qili Shen
Kairui Feng
Ang He
Xuanhong Chen
Dake Zhang
Venue: ICLR 2026 Poster (2026)
Keywords: AI for Science
Reinforcement Learning
PDE Control
Surrogate Modeling
Primary area: reinforcement learning

Abstract: Recent work increasingly adopts AI for Science (AI4S) models to replace expensive PDE solvers as simulation environments for reinforcement learning (RL), enabling faster training in complex physical control tasks. However, using approximate simulators introduces modeling errors that affect the learned policy. In this paper, we introduce a unified $\varepsilon$-$N$ framework that quantifies the minimal computational cost $N^*(\varepsilon)$ required for an AI4S model to ensure that tabular RL can estimate the value function with unbiasedness, with probability at least $1 - \delta$. This characterization allows us to connect surrogate accuracy, grid resolution, and RL policy quality under a shared probabilistic language.  We analyze how the discretization level $K$ of AI4S and RL space governs both PDE surrogate error and RL lattice approximation error, and we employ spectral theory and Sobolev estimates to derive optimal grid strategies that minimize total cost while preserving learning fidelity. Our theory reveals that different systems, such as ODE- and PDE-governed environments, require different allocations of effort between physical simulation and RL optimization. Overall, our framework offers a principled foundation for designing efficient, scalable, and cost-aware AI4S-RL systems with provable learning guarantees.

Published as a conference paper at ICLR 2026
ON THE COMPUTATIONAL LIMITS OF AI4S-RL A UNIFIED E-N ANALYSIS
Qili Shen Kairui "oogl..-Ann He,, Xuanhong Chen Dake Zhang
TTngji University, Shanghai, China
Shanghai Innovation Institute, Shanghai, China (2311788 kelviiiinn kelvinfkr' @tonggii edu,
*Shanghai Jiao Tong University. Shanghai, China angh@avidia. com, (chenxuanhong, dk. zhang] edu. en
ABSTRACT
Recent work increasingly adopts AI for Science (AI4S) models to replace expensive PDE solvers as simulation environments for reinforcement learning (RL., enabling faster training in complex physical control tasks. However, using approximate sim- ulators introduces modeling errors that affect the learned policy. In this paper, we introduce unified e-N framework that quantifies the minimal computational cost (*) required for an AI4S model to ensure that tabular RL can estimate the value function with unbiasedness, with probability at least - 6. This characterization allows us to connect surrogate accuracy, grid resolution, and RL policy quality under shared probabilistic language. We analyze how the discretization level K of A14S and RL space governs both PDE surrogate error and RI. lattice approx inptionmation error, and we employ spectral theory and Sobolev estimates derive optimal grid strategies that minimize total cost while preserving learning fidelity. Our theory reveals that different systems such as ODE- and PDE-governed environments require different allocations of effort between physical simulation and RL optimiza- tion, which is consistent with the empirical results. Overall, our framework offers a principled foundation for designing efficient, scalable, and cost-aware A14S-RL systems with provable learning guarantees.
INTRODUCTION
Reinforcement learning (RL) for PDE-constrained control promises transformative impact in domains such as fusion energy (Degrave et al.) 2022) and climate modeling (Feng et al. (025). Yet fundamental obstacle remains: high-fidelity PDE simulations are far too costly to support the millions of interactions demanded by RL. A natural remedy is to replace PDEs with Al surrogates, but this introduces critical and largely unexplored question: how should the resolution of the surrogate be coordinated with the RL agent's discretization as to minimize total computational cost while preserving policy accuracy?
Classical RL theory assumes stochastic transition noise that vanishes under repeated sampling (Sutton al. 1998 Azar et al. (017),, and robust RL extends this to worst-case distributional uncertainty (Der man et al. 2021 Agarwal & Zhang 2022). In contrast, AI4S surrogates incur deterministic errors stemming from spatial discretization, temporal integration, and boundary approximation. Such errors do not average with more data and can systematically bias policy learning. In a cart-pole system, if the surrogate discretization does not adequately capture the upright equilibrium, the learned controller will consistently overshoot; regardless of how fine the RL action resolution is, policy accuracy is ultimately constrained by the surrogate's 's resolution. Intuitively, minimizing regret under a fixed computational budget requires balancing the grid resolution of the AI4S surrogate with the decision resolution of the RL system.
"Equal contribution. Corresponding authors.

---

Published as a conference paper at ICLR 2026
To address this trade-off, we introduce framework that explicitly models discretization-inuued bias. By sampling multiple trajectories from perturbed initial conditions within observation uncertainty bounds, recast deterministic surrogate eπrors statistical inference problem. This structured sampling bridges numerical analysis in A14S with RI. learning theory, enabling principled resolution coordination between surrogates and agents. We focus on tabular RI. to establish necessary condition for learnability, providing a resolution lower bound that guarantees feasibility and guides hyperparameter design for deep RL in AI4S systems.
Contributions. We introduce the e-N framework that quantifies the minimal number of surrogate calls N required to achieve e-accurate value estimation with confidence 1- -. assuming optimal resolution allocation between AI4S simulation and RL exploration. Our analysis reveals that this optimization problem exhibits surprising complexity: optimal resolution ratios follow non-obvious fractional scaling laws (e.g.. At, = (At..) )2/3 for tokamak control, where At, is the discrete timestep in the RL domain and a is that in the A14S/World domain), have exponential sensitivity to system dynamics, and vary significantly across different physical systems. We provide:
1. Resolution coupling theory: Our p-K analysis shows that the forward prediction error p scales as OL/HH, + 1 in the high-resolution regime, where H,, is the RL's temporal resolution (inverse decision frequency), K is the ratio of physical grid cells to RL states, and d is the spatial dimension of PDE system. This result quantifies how the resolutions of both the AI4S surrogate and the RL agent jointly determine the learning quality.
We derive system-specific expressions for the optimal res- olution ratio K* that minimize computational cost. For instance, K = (7/4 - exp(1/A1)) for tokamak plasma control and K = 2 exp(1/A1) for cart-pole systems, where a denotes the dominant modal growth rate of the PDE dynamics.
Closed-form optimal allocations:
3. Empirical validation: Experiments confirm that both tabular RL and DQN exhibit the predicted non-monotonic sensitivity to resolution, with off-optimal parameter choices requiring computational cost scaling as approximately N relative to the optimal N to reach the same accuracy. Moreover, there exists an optimal grid ratio K and an optimal balance between H, and H.., and these values are found to be relatively consistent across different learning algorithms.
Related Literature. Our work lies the intersection of reinforcement learning theory and scientific computing. On the RL side, the tabular literature has established Probably Approximately Correct (PAC) bounds that characterize sample complexity (Sutton et al., 1998 Azar al.) 2017 Dann et al.
2017). Extensions to function approximation further reveal how this complexity depends structural properties such as Bellman rank and eluder dimension (Jin et al. 2020]] Wang et al. 2020 Jin et al., 2021 A key limitation of these analyses is their reliance on stochastic transition dynamics, where errors are assumed to average out. This assumption breaks down in scientific computing surrogates, where deterministic discretization errors introduce systematic, resolution-dependee
biases. While robust RL frameworks (Derman et al. 2022 AAgaraal & Zhang, [022] Shi et al.) 2023 Kwon et al., 2021 address model uncertainty through worst-case optimization, they primarily target irreducible stochastic noise. In contrast, our focus is on structured numerical biases, which we exploit to derive principled coallution-sampll trade-offs.
On the scientific computing side, recent advances in neural solvers have accelerated PDE simulations Lu et al., [002 Kovachki et al. 2023)
by orders of magnitude. Neural operators (Li et al.) 2020
and physics-informed neural networks (Raissi et al. 2019 Karniadakis et al. 2021) have enabled
breakthroughs weather prediction (Bi et al. 2023), turbulence modeling (Kochkov et al. 202 )), and
plasma control (Degrave et al. 2022 Although operator learning admits error bounds (Kovachki et al. 2023), prior work rarely investigates how surrogate resolution affects downstream RL performance. Classical PDE control methods rely on adjoints and Pontryagin's 's principle (Pontryagin 2018), 3,, and recent applications of RL to PDEs (Farahmand et al. 2017 Han et al. 2018) typically assume access to true dynamics or ignore surrogate emror..
2

---

Published as a conference paper at ICLR 2026
2 PRELIMINARIES
Notation. We study a hybrid system where an RL agent interacts with an A14S surrogate approxi- mating a PDE-govemed environment. comprehensive notation table provided in Appendix ^ with key notations:
Subscripts r and w denote RL and physical world parameters respectively.
K = Ar,/ A is the ratio of spatial resolutions between the RL and A14S domains, where Ar, the spatial grid size in the RL domain and Arw w that in the AI4S domain.
Hr and H. denote the temporal discretization (time-step size) of the RL agent and the A14S surrogate, respectively.
E, 5 are the target accuracy and confidence levels.
AI4S Surrogate Errors. A14S models approximate PDE solvers through neural operators. For PDE solution operator G - approximated by surrogate Go, the error decomposes as:
[ - Go(All S C, + C2X
(2.1)
where H is temporal discretization, X is spatial resolution ratio, and s, d are parameters determined by the dimension and property of the underlying PDE solver.
Probably Approximately Correct (PAC) Learning in Tabular RI. We define episodic Markov Decision Process (MDP) as M = (S, A, P, R. H), where S is a finite state space with [S] = S, A finite action space with |A| A, P is an unknown transition kernel, R A [0, 1] the reward function, and H is the episode horizon. For clarity, we omit the subscript r here; in later sections, we will use distinguish RL parameters from those of the underlying physical system.
We call that an algorithm is PAC-MDP if for any ō E (0, 1), with probability at least 1 8.
outputs policy satisfying Vi (sp) V (sn) S after most poly(S. A, H. 1/E, log (1/8)
episodes. A classical result in tabular RI. settings is that the UCB-VI algorithm Azar et al. (2017)
achieves PAC-MDP guarantee with sample complexity Nepisodes Ö |AAf" log } logarithmic factors.
.where hides
EMPIRICAL MOTIVATION: RESOLUTION COORDINATION IN AI4S-RL
Before developing our thcoretical framework, we first empirically demonstrate that an optimal discretization scale exists in AI4S-RL systems. This balance between surrogate and agent resolutions follows scaling relation across grids.
We utilize the simulated Cart-Pole environment, which is tractable for exhaustive search yet represen- tative of continuous dynamical systems. We establish two distinct experimental settings to bridge theory and practice. First, we employ Tabular Value Iteration the exact physics engine a theoretical oracle. Second, to strictly emulate the AI4S workflow where learned dynamics replace expensive solvers, we train Deep Q-Networks (DQN) within surrogate environment governed by pre-tra are-trained neural network. Across both settings, we vary the temporal resolution ratio log H, and the spatial resolution ratio K, measuring the sample complexity required achieve value function error below 1 %.. Full experimental details are provided in Appendix C.5.3
as
Figures D and 2 show that sample complexity does not improve monotonically with resolution.
Instead, both methods exhibit similar optima: tabular RL at (K = 1.5, log H, 1/3), and Q-learning at (K 2.0, log H, = 1/2). Off-optimal parameter choices require computational cost scaling as approximately 1.6 relative to the optimal N order to reach the same accuracy.
These findings highlight that conventional hyperparameter search strategies (e.g., grid or binary search) are ineffective in this setting: resolution parameters span multiple orders of magnitude, making brute- force exploration prohibitively expensive in high-dimensiioal PDE systems. This motivates the development of a theoretical framework that predicts optimal resolution trade-offs from system properties, rather than relying on empirical tuning.

---

Published as a conference paper at ICLR 2026
Sample Complesity Heatmep
Q-leaminn Sample Complesity Heatmmg
6:11
4:1
5.20
244
4.5
4.5
423
125
2/5
2.8
472
241
31
15
4.10
4.8
415
634
1/3
1:50
2.5
2.68
211
10
1/e
1/9
10
4.25
4.3
422
436
2.05
4.5
4/0
4.2
466
5,50
2.5
1.00
361
27
CVetu
Figure 1: Sample complexity heatmap for the Cart-Pole system under tabular RL, showing non- monotonic dependence on resolution parame- ters. The red box indicates optimal configura- tions. Color intensity corresponds (NV as shown in the colorbar.
Figure 2: Sample complexity heatmap for the Cart-Pole system under Q-Learning, revealing algorithn-specific optimal configurations. Theo- retical analysis is provided in Section 4.4
4 ERROR PROPAGATION AND RESOLUTION COORDINATION IN A14S-RL SYSTEMS
In this section, we develop theoretical framework that predicts optimal resolution trade-offs from system properties. In Section 4.1 we show that the discretization the AI4S surrogate itself determines whether the RL agent converge PAC guarantees. In Section 4.2 analyze how the interaction between AI4S and RL discretization affects the convergence rate of the state transition matrix, using magnetohydrodynani (MHD) tokamak as a representative example. Building on these error models, Section 4.3 investigates the trade-off between the computational costs of surrogate construction and RL training. We derive the optimal resource allocation strategy- specifically the discretization resolutions that minimizes the total cost required to achieve PAC guarantee. Finally, in Section 4.4] we compute the optimal discretization and computational cost in four real-world systems to validate the framework.
4.1 SPECTRAL ANALYSIS OF ERROR AMPLIFICATION IN
A14S SYSTEMS
In AI4S systems governed by physical laws (e.g., PDEs), predictions are obtained by evolving the current state yo through a deterministic solver f under action a, yielding y1 = f a). However, yo
is only observed up grid size Ay, so the predicted next state lies in the range [y] y/ y/ where the bounds are defined by yi mii = / f (3/ + Ayo. a) and ymax = max (yy f Ayo, a). Here, Atjo the measurement uncertainty in practice.
If the interval [yjmin y]ma]] spans more than two full grid cells of PDE systems, the next state cannot be uniquely identified, and PAC guarantees are unattainable. In contrast, the interval covers only single full grid cell (with any additional coverage being partial), repeated sampling of the initial state allows us to estimate the next PDE state with arbitrarily high confidence. The convergence rate of this estimation is summarized in Theorem
Theorem 1 (Sample Complexity for & d-Confidence Classification). Consider repeated forward predictions from perturbed initial states, where the predicted frequency is the empirical probability of the next state falling into a given grid cell. Let p denote the predicted frequency with which the true next state falls into the correct grid cell, and let q = Pmax denote the maximum predicted frequency among all competing cells. To resolve the correct cell with confidence at least 5, the number of
forward predictions required is bounded by n
log(1/8)
min: Agy(a) Panaa
4

---

Published as a conference paper at ICLR 2026
Proof in Appendix C.I.7
Remark 1 (Time-Step Constraint for State Separability). Consider the solution operator f<( of nonlinear PDE. For observation perturbations y with ||d|| S Ay, the linearized perturbation dynamics admit a modal decomposition nyt) ~ Ex fie where {Ak} ) are the eigenvalues characterizing modal growth rates. Classification-cell separability requires that the dominant perturbation growth remains bounded, ie., e I 1, where A == max (Re(Ae)) This implies a time-step constraint At S 1/21.
Proof in Appendix C.1.2
From Theorem when the two largest cell frequencies Ay6) and piilx are close, where (Ay() Pmaxx )
(5) )² 0 and o (Aghd)
log(1/8)
oc, PAC learning becomes unattainable. Remark
shows that smaller At enlarges the frequency gap (Ayli) Pmax),, while if At 1/A1, perturbations exceed Ay and cell frequencies equalize, destroying state separability. Hence the constraint A/ S 1/A is necessary for PAC guarantees. We note that defining the global maximum is primarily for theoretical parsimony; utilizing a step-specific local growth rate in practice would not fundamentally alter the complexity analysis.
4.2 ERROR COUPLING BETWEEN RL AND PDE SPACES IN AI4S SYSTEMS
showed that the intrinsic dynamics and temporal discretization of the physical system determine whether PAC guarantees are achievable, and computed the sample complexity to reach PAC when possible.
In Section 4.
In this section, we extend the analysis to mismatched discretizations between the RL and PDE spaces. Each RI. action a, must be lifted to the PDE space for forward evolution, and the resulting physical state projected back for RI. These projections introduce errors, which may be amplified near PDE boundaries by nonlinear operators. Since RL actions often operate directly on boundaries (c.g.. controlling droplet interfaces in tokamaks to avoid wall contact, shaping fluid boundaries around an aircraft wing to increase lift), boundary accuracy crucial. In MHD systems with droplet interfaces, the trace theorem implies that poor boundary approximation further enlarges these errors. To clarify the relationship between the two spaces, we define the discretization tuples (H,, S,, A,) for the RL agent and (Hw, wr, Sw., Au:) for the PDE system, where the components correspond temporal, state, and action (or control parameter) granularities, respectively.
a. RL Space
b. Tokamak Device
c. PDE Space
E. MHD governing equations
Isoflux
Line
Dr (pv)
Vessel
(Continuons equation)
Plasma
Vy Vp pVv
Faxtss
position
(Momentum equation)
Strike
8B at
points
V (v B) + "
(Magnetic field control equation)
V 0
(Deformttion the droplet's interface)
o% = Trattes (Boundary conditions)
d. Surface Interpolation
Interpolation Error
Figure 3: Illustration of Error Propagation in a Tokamak Machine RL controlling System.
MHD System As representative case, this section we focus on a magnetohydrodynanme (MHD) tokamak control task, where an RI. agent manipulates boundary magnetic coils to suspend plasma
5

---

Published as a conference paper at ICLR 2026
droplet at the chamber center without wall contact Degrave al. (2022); Ding et al. (2024). We assume observations of the MHD fluid and magnetic field on a predefined spatial lattice, and conduct numerical validation based realistic tokamak configurations.
Figure 3 illustrates the governing PDEs, interface tracking methods, and control structure involved in the RL.bbased regulation of the droplet. The reward signal for reinforcement learning is defined the minimum distance between the plasma droplet and the tokamak boundary, and the control actions correspond to adjustments of the magnetic field via coil currents. These equations define the error sources and discretization structure that form the basis for our multi-resolution analysis.
We assume continuous control over the magnetic field via time-dependent coil currents (1), enabling real-time adjustments to maintain droplet suspension without wall contact. Reinforcement learning algorithms adapt these control imputs based on observations of the droplet boundary and plasma dynamics. Key physical variables variable--including the magnetic field B, current density J, density p. velocity field v, and interface function Φ are continuously monitored and regulated to optimize positioning and stability.
Interface Error Propagation via described above, RL actions do not directly affect the plasma interface but propagate indirectly from boundary coils through multiple physical fields. This indirect coupling introduces compounding discretization errors. We first analyze ainle-equation timestep errors, then demonstrate their coupling across the full system.
Coupling In the tokamak control setting
The control-to-interface coupling follows a hierarchical propagation path:
RL action d,
projectinn
Adw
houndary,
B an
curl
Lorenzz
advection
(4.1)
where each arrow represents a potential error amplification point. The magnetic field B controlled at the boundary generates current density J V x B, which produces Lorentz forces that drive the velocity field v, ultimately advecting the interface level set Φ.
At each stage, discretization emror accumulate. The RL action a, must be mapped to physical control parameters aw. introducing projection error (Aa, - Aaw). Magnetic field boundary conditions suffer from reduced regularity, yielding [B B, (aa) = O(I Axp.aal (//p/all//)) by the trace theorem.
These control errors then propagate nonlinearly into the velocity field as Av = O affecting the entire flow dynamics.
06 The droplet interface evolves according to the level set equation DE + Euler time discretization and central spatial differences, the local truncation error per PDE timestep
VΦ O. Under forward
becomes (kk) ,tt,,,)) ) + O(I ((,)x,
However, incorporating the velocity perturbations from boundary control and the reduced regularity interfaces, the complete per-step interface error is:
A (k)
((Stw) temporal
O(1AAx piml)) mlll)
O(I (XIAxx, 1/2) + O(Da |1Axp, 2/2)
spatial interior
control
boundary
coupled
(4.2)
In the following paragraph, we will couple this error with contributions from other governing equations to enable a comprehensive analysis.
Refined Total Error Decomposition with Time-Scale Separation By extending the interface error strueture discussed previously, generalize the analysis to broader class of PDE-based surrogate environments (details in Supplementary Information). Over a single RL step, the total prediction error can be decomposed as Awai C2 || + C.2
Os Afr (I Ax,, + A 1/2 || 1/2 Aa 1/2 AI papppall//))) /12). where the first four terms represent RL-side errors (interior space, boundary space, action space, and temporal discretization), and the last term captures PDE surrogate errors amplified over the RL-to-physies time-scale ratio At,/ /t
To ensure that the RL transition kernel remains distinguishable under observation uncertainty Ay, we impose the condition Awat O(Ay). This leads to the following resolution matching constraints

---

Published as a conference paper at ICLR 2026
between RL and AI4S components:
| || ~
A int ~ Aat/2 Aa Sl. ~ Ay,
(4.3)
|A
~
Axp.tal ~ Av?, Atw < Ay.
This decomposition makes explicit how discretization scales in the RL and AT4S pipelines must be jointly selected to ensure unbiased value estimation. In particular, the term /t, At/ highlights how fine-scale physical errors can be amplified over coarser RL horizons, motivating resolution-aware control design.
In Section 4.1 we analyzed how to approximate the RL transition kernel under observation noise Ay. Here, we incorporate numerical emrors from the AI4S prediction and define the misclassification rate after one RI-to-AI4S projection and back. Concretely, let denote the probability that the true next physical state given observed state falls into the correct grid cell after (i) projecting the observed PDE state into the RL grid, (ii) selecting the truncated optimal RL action, and (iii) mapping this action back the PDE. Then P 1-p. which quantifies the probability that one RL-PDE interaction fails to preserve grid-level consistency.
From Eq. 4.3 total numerical error Amnal ~ &t ~ Ey and letting the RL temporal resolution be H,, T, (we also use this further on unit for calculating total computational cost), the total error rate is:
Ay p=l- \/Ay/ H Ay Atntal
(4.4)
X1/H, + C1/K3
Hence, in the high-resolution limit (when H, and K3 are large) p= 0= finer RI. temporal resolution and A14S spatial refinement both reduce the forward projection error rate at inverse polynomial rates.
showing that
From more general analysis on PDE scales, we could get:
Theorem 2 (p-K Analysis for d-Dimensional Systems). For a d-dimensional PDE system, the forward projection error rate is: p = A1HH+++C.//k where the numerical error scales as
A Amial K-K Byy from spatial discretization. In the high-resolution limit: p p= 0= K
Proof in Appendix C.I.4
4.3 OPTIMAL COMPUTATIONAL COST ALLOCATION BETWEEN RL AND A14S
To translate these PAC-possible A14S-RL analysis into practical system design guidance, this section establishes computational resource allocation framework for A14S-RL systems. AI4S models are generally trained under fixed computational budgets, we aim to derive how their discretization parameters dentted as Hw. Sw., Aw-can can be aligned with the RL-side parameters H,, Sp, Ar to ensure computational and statistical consistency.
can
On the physical side, computational cost scales as H., - For tabular RL, classical UCB- VI methods require sample complexity of ((HA,,A,) for value function convergence without any assumptions on transition smoothness. Based on our transition identifiability analysis, an improved bound of OH,SS,
A...
miny (Ay()) pla). 5²) )²)
log(1/5)
achievable when the transition matrix is
known up to statistical confidence 8, reducing the RL problem to dynamic programming. From our p-K analysis, we approximate pinax & pAy
т,, - Ay, leading refined sample complexity of O (HASA log(1/8)) which we adopt in the remainder of this section.
(f)
We then align RL and AI4S computational costs via the relation S A., = H A,. From the interface and baidd error analysis in Section 4.2 for MHD system, and the temporal resolution bound in Section 4.1 we discretize the physical space in three spatial dimensions (æ, 2), where:
A, OAA.) (z²)), S, O(S.) ( ), H, O(H2/3), H.. O S x This constraint also respects CFL stability conditions for fluid systems.
7

---

Published as a conference paper at ICLR 2026
By substituting these scaling relations into cost-balancing equation, express the total computational cost as a function of the resolution ratio K, leading to the following theorem.
Scaling). For a physical system subject to projection error with state space scaling S, K and action space scaling A, ~ K the com- putational balance condition H. A, H. Sw. ensures equivalence between the RL surrogate resolution and the underlying PDE world resolution. Minimizing the overall computational costs,
Theorem 3 (Optimal Resolution with gystem-Dependee
(103(1/0)
(log(1/0) E2
Cost( K) = H2K
(4.5)
H,
yields optimal resolution ratio between RL and AI4S space:
a ++ 2d
+
+
2d
(a + (3 -H K
=
exp
(4.6)
a+
a
dYp
when H, 2 11
Proof in Appendix C.I.5
Remark (Impact of Actuation Topology). The scaling exponents derive from the Trace Theorem based on the control location. Boundary-actuated systems (e.g., tokamak) require quadratic refine- ment due boundary regularity loss, yielding systems (e.g., heat sequencing) follow follow standard volumetric scaling with a = d and 3 = da- where d denotes the physical spatial dimension and da the action dimension. This distinction dictates the optimal resolution K* and the cost disparities shown in Table
2d and 3 2da- Conversely, interior-actuated
4.4 THEORETICAL ANALYSIS ACROSS PHYSICAL SYSTEMS
Building on our theoretical framework, now evaluate how surrogate discretization, spectral response, and RL precision jointly affect sample complexity across four representative A14S RL systems: (1) Tokamak plasma control, (2) Turbulent airfoil regulation, (3) Teppanyaki heat sequencing. and (4) Cart-pole stabilization (Details in Supplementary). Figure iillssrrtes their discretized state spaces.
Optimal vs
Cart-Pole
Intorin Mach
Argle
Cart Fale Syown
Marseic
Teppanyaki
Tokamak
Figure 4: State space discretization for four A14S-RL systems.
Figure 5: Log-scale plot of K* vs. modal growth X
We adopted an e-N perspective: to ensure value function estimation within error & and confidence 1 1- 5, what is the minimal computational cost N, given surrogate discretization and physical uncertainty? Based on our analysis, we can now systematically achieve optimal resolution allocation under (c, 8) accuracy constraints through unified framework that integrates normalization, identifiability analysis, and closed-for cost expressions ppphess minimizing the total computational cost of AI4S-RL systems while ensuring reliable policy estimation.
Table I summarizes scaling relations and optimal resolution ratios. Systems with strong boundary observability (e.g., tokamak, airfoil) require quadratic scaling in state-action discretization relative
8

---

Published as a conference paper at ICLR 2026
to mesh resolution, while low-dimensional systems tolerate aggressive upsampling without cost explosion. Figure [ shows how the optimal K' saturates with increasing λ1, especially in high- dimensional PDEs. In contrast, ODE-based systems benefit more directly from higher surrogate
resolution.
Our theoretical framework provides insights into the empirical results from Figures and For tabular RL, the observed optimum at (K = 1.5, logn H, /3) requires 103.66 2s 4500 samples, while the worst configuration at (K = 2.5, logH H, = 2) needs 10 as 46000 samples-a 10-fold increase. This translates to computational cost scaling approximately as N when deviating from optimal parameters, consistent with our claim in the abstract. Q-learning achieves its minimum at (K 2.0, log H, = 1/2) with only 102-65 = 450 samples, demonstrating that function approximation can effectively exploit state similarity to reduce sample complexity by an order of magnitude.
Table Optimal Resolution Matching and Computation Cost Across Systems under error of value function estimation and
5 confidence level
System
Resolution Scaling
Computational Cost
Optimal K*
Tokamak Control
S. OSS,) ((2
A. = O(A..) = (
Н. HD O(sw) < <//A1
H²K*
22
(1
(1/(1/6)
H.
K3 )z
i
1/3
exp
311
A, O(A..) =
K² 6
log(1/6)
1/2
Airfoil Control
S. O(S,) = 0(x4))
H, H O(zw) 1/A1
2²
xCz )2 (1 TT.
exp
2
Teppanyaki Plate
S. ((S.,) = 0(x
A, O(A.) ((sw)
H²K
log(1/8)
7
1/2
1
H. = H (((u)) < <1/A1
(1 TI. R2
3
2A1
A. O(A..) O(Iw)
H²K²
log(1/0)
Cart-Pole System
S. O(S..) = -(ww)
2. exp
H. /3 O(ww) < 1//11
c²
(1
A1
H.
These results suggest that systems with high spectral amplification yield diminishing returns from A14S refinement unless RL granularity is co-optimized. Our s-N framework thus provides principled tool for balancing simulation fidelity and RL efficiency under finite compute budgets.
4.5 EMPIRICAL VALIDATION ON PDE- IDE-GOVERNED CONTROL: TEPPANYAKI HEAT SEQUENCING
While the Cart-Pole experiments Section 3 successfully demonstrate resolution trade-offs, they rely low- dimeenenioonl ODE system with learned dynamics. To rigorously validate theoretical predictions in sigh-dimensionall PDE environments with continuous state-action spaces, we conduct comprehensive study a two-dimensional control problem: multi-item heat sequencing teppanyaki cooking surface. This environment is governed by the canonical two-dimensional heat diffusion equation ØT /Ot a 2 + Q t), representing a broad class of parabolic PDEs encoun- tered in A14S applications (thermal management, plasma diffusion, etc.). The control task requires simultaneously cooking three food items at distinct temperatures using two heat sources, posing challenges for policy learning under discretization constraints. We train Proximal Policy Optimization (PPO) agents across Cartesian product of spatial resolution ratios K and temporal resolution param- eters log, (hr,), yielding 25 experimental configurations as detailed in Figure 6 All experiments run on identical hardware (Intel Xeon Gold 6530, NVIDIA RTX 4090). Complete environment specifications, PPO hyperparameters, and convergence criteria provided in Appendix E
Figure 6 illustrates the comprehensive experimental results on the Teppanyaki thermal control task. Panels (a) and (b) present the environment setup: the agent regulates distributed heat sources to manage non-stationary diffusion dynamics for precise multi-object temperature tracking. Panel (c) displays E-N scaling curves for all 23 qualified configurations, where color denotes temporal resolution and intensity indicates spatial resolution K. Panels (d), (c), and (f) show CPU, GPU, and total computational cost heatmaps across the (K, y) configuration space. Configurations at K 8.0

---

Published as a conference paper at ICLR 2026
3. Optimal Control results
b. Teppanyak Environment
e-N Scating Curves
87 =av/'' 4/000
Loo Sule
guffft 2000
Salid regs Rooa
Trilking Sago
3-
155
12
25
22
27
201
351
150
15
18
22
27
381
21
-17
115
201
17
15
122
125
3
220
12
12
13
d. GPU Computational Cost
CPU Computationl Coat
(.Toaal Computationl Coat
Figure 6: Experimental analysis of the Teppanyaki PDE control task. (a) The analytical optimal control solution serves as a performance oracle, illustrating the precise temperature tracking required. (b) Visualization of the simulation environment and the governing MHD/heat diffusion equations. (c-e) e-N scaling curves (top) with CPU/GPU computational cost decomposition heatmaps (bottom), demonstrating the non- non-monotonic resolution-effciieng trade-off.
(marked "N/A") fail to converge, validating that excessive coarsening violates learnability conditions. The optimal configuration (K* = 6, log (hz,) 1/3) achieves minimum total cost, confirming that Theorem provides effective guidance for deep RL hyperparameter selection.
5 CONCLUSION AND LIMITATIONS
This work presents unified theoretical framework for analyzing RL in surrogate-based physical environments governed by PDEs. By integrating spectral discrctization theory with PAC sample complexity bounds from tabular RL, we establish how surrogate resolution, physical dynamics, and RL precision jointly influence the number of episodes required for accurate policy learning.
We propose the E-N framework characterize the minimal RL-AI4S interaction cost under fixed accuracy and confidence requirements. Through analysis of surrogate-induced error propagations, and optimal resolution trade-off, we derive system-specific cost scaling laws and identify spectral regimes where surrogate refinement provides diminishing returns. Empirical validation across four representative AI4S-RL systems confirms that spectral growth rate, boundary observability, and dimensionality all play critical roles in determining learning efficiency.
Several assumptions in our analysis warrant further investigation. First, our analysis focuses on tabular RL for theoretical tractability and to isolate physical discretization errors, establishing necessary condition for learnability that applies to any algorithm. While our experiments with DQN and PPO confirm that these resolution trade-offs govern the primary performance trends even in deep RI., rigorously extending our mathematical bounds to function approximation settings remains an open direction. Second, treat surrogate models fixed and externally trained; adaptive refinement or active error correction within the RL loop could further improve performance. Finally, framework focuses deterministic discretization errors. While real-world neural surrogates also introduce approximation and generalization errors, our analysis establishes a theoretical baseline. isolating the resolution trade-off as a fundamental constraint.
Overall, this work provides a principled foundation for evaluating and co-designing A14S surrogates and RL policies in resource-constrained scientific applications.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Alekh Agarwal and Tong Zhang. Model-based rl with optimistic posterior sampling: Structural conditions and sample complexity. Advances Neural Information Processing Systems, 35: 35284-35297, 2022.
Mohammad Gheshlaghi Azar, lan Osband, and Rémi Munos. Minimax regret bounds for reinforce- ment learning. In International conference on machine learning, pp. 263-273377 PMLR, 2017.
Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate medium- range global weather forecasting with 3d neural networks. Nature, 619(7970):533-533 2023.
Christoph Dann, Tor Lattimore, and Emma Brunskill. Unifying pac and regret: Uniform pac bounds for episodic reinforcement learning. Advances in Neural Information Processing Systems, 30, 2017.
Jonas Degrave, Federico Felici, Jonas Buchli, et al. Magnetic control of tokamak plasmas through deep reinforcement learning. Nature, 602(7895):4144499) 2022.
Esther Derman, Matthieu Geist, and Shie Mannor. Twice regularized mdps and the equivalence between robustness and regularization. Advances in Neural Information Processing Systems, 34: 22274-22287, 2021.
Shiping Ding, AM Garofalo, HQ Wang, DB Weisberg, ZY Li, X Jian, David Eldon, BS Victor, Alessandro Marinoni, QM Hu, et al. A high-density and nigh-confinement tokamak plasma regime for fusion energy. Nature, 629(8012):555-560) 2024.
Weinan E. Principles of Multiscale Modeling. Cambridge University Press, 2011.
Lawrence C Evans. Purtial differential equations, volume 19, American mathematical society, 2022.
Amir-massoud Farahmand. Salch Nabi, and Daniel N Nikovski. Deep reinforcement learning for partial differential equation control. In 2017 American Control Conference (ACC), pp. 3120-3127. IEEE, 2017.
Kairui Feng, Ning Lin, Robert E Kopp, Siyuan Xian, and Michael Oppenheimer. Reinforcement learning-based adaptive strategies for climate change adaptation: An application for coastal flood risk management. Proceedings of the National Academy of Sciences, 122(12):c2402822122 2025.
Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115((4))8505-88 -510008500,,,, 2018.
Thomas Y Hou. Numerical approximations to multiscale solutions in partial differential equations.
Frontiers in Numerical Analysis:: Durham 2002, pp. 241-301, 2003.
Chi Jin, Zhuoran Yang, Zhaoran Wang, and Michael Jordan. Provably efficient reinforcement learning with linear function approximation. In Conference on learning theory, pp. 2137-2143, PMLR, 2020.
Chi Jin, Qinghua Liu, and Sobhan Miryoosefi. Bellman eluder dimension: New rich classes of rl problems, and sample-efficient algorithms. Advances in neural information processing systems, 1113406-13418. 2021.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6))422-440. 2021.
Dmitrii Kochkov, Jamie Smith, et al. Machine leamning-accelerateed computational fluid dynamics. Proceedings of the National Academy of Sciences, (1802)):e210178411 2021.
Felix Kochler, Simon Niedermayr, Nils Thuerey, al. Apebench: A benchmark for autoregressive neural emulators of pdes. Advances in Neural Information Processing Systems, 37:120252-122 120310, 2024.
11

---

Published as a conference paper at ICLR 2026
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadeneshelii Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1-97) 2023.
Jeongyeol Kwon, Yonathan Efroni, Constantine Caramanis, and Shie Mannor. RI for latent mdps: Regret guarantees and a lower bound. Advances in Neural Information Processing Systems, 34: 24523-24534, 2021.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.0888 00800, 2020.
Jacques-Louis Lions. Quelques méthodes résolution des problèmes aux limites non linéaires. Dunod; Gauthier- Villaas,, 1969.
Lu Lu, Pengzhan Jin, Guofei Pang. Zhen Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218 229, 2021.
Andrew Majda, Andrea I. Bertozzi, and A Ogawa. Vorticity and incompressible flow. cambridge texts in applied mathematics. Appl. Mech. Rev, 55(4):B77-B78, 2002.
Savinay Nagendra, Nikhil Podila, Rashmi Ugarakhod, and Koshy George. Comparison of reinforce- ment learning algorithms applied to the cart-pole problem. In 2017 international conference on advances in computing, communications and informatics (ICACCI), pp. 26-32. IERE, 2017.
Lev Semenovich Pontryagin. Mathematical theory of optimal processes. Routledge, 2018
Koldo Portal-Porras, Unai Fernandez-Gamiz, Ekaitz Zulueta, Alejandro Ballesteros-Coll, and Asier Zulueta. Cnn-based flow control device modelling aerodynamic airfoils. Scientific Reports, 12 (1):8205, 2022.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Laixi Shi, Gen Li, Yuting Wei, Yuxin Chen, Matthicu Geist, and Yuejic Chi. The curious price of distributional robustness reinforeement learning with a generative model. Advances Neural Information Processing Systems, 36:79903-7999 -799177 2023.
Richard S Sutton, Andrew G Barto, et al. Reinforcement learning: An introduction, volume 1. MIT press Cambridge, 1998.
Ruosong Wang, Russ Salakhutdinov, and Lin Yang. Reinforcement learning with general value function approximation: Provably efficient approach via bounded eluder dimension. Advances in Neural Information Processing Systems, 33:6123-6135, 2020.
12

---

Published as a conference paper at ICLR 2026
A TABLE OF NOTATIONS
The following table provides reference for the key symbols and notations used throughout this paper.
Symbol
Description
General Reinforcement Learning (RL) Parameters
S,A
The finite sets representing the state space and action space.
H
The horizon, representing the length of an episode.
[AA
The size of the state and action spaces,
The desired error tolerance for value function estimation. The confidence parameter.
(T(s)
The value function of a policy и at state s.
V
The optimal value function.
The optimal policy and the learned policy, respectively.
P.P
The true and the empirically estimated state transition kernels.
The sample complexity.
N
The overall computational cost is treated as a function N our e-N analysis.
AI4S-RL Framework Parameters
Sr, Ar, Hr
The scale of the state space, action space, and horizon in the RL domain. The scale of the state space, action space, and horizon in the AI4S domain. The discrete timesteps in the RL and AI4S/World domains, respectively. The spatial grid resolution the RL and A14S/PDE domains. The Resolution Ratio, defined K A Ae
Sw: Am, Hy
At,, Atw
Ary, Atw
K
K
The optimal resolution ratio that minimizes computational cost. An initial state and the measurement uncertainty in practice.
Vo: Ayo
Ay
The grid cell size.
Atotal
The total numerical error accumulated in RL-AI4S interaction step.
The forward projection eтror rate in the AI4S RL system.
PDE and Spectral Analysis Parameters
A1
The dominant (largest real part) eigenvalue of the linearized PDE operator. The Fréchet-linearized operator of the PDE system around state yo- The k-tth orthonormal basis function.
C(yo)
ők
7x
The modal gain factor for the k-th mode at time t.
(t)
The propagated mode shape for the k-th mode at time t.
Mathematical Spaces and Operators
(
The Sobolev space of order over a domain Ω.
IP(R)
The Lebesgue space of p-integrable functions over domain Ω.
V,A
The gradient and Laplacian operators, respectively.
-lv
The norm in vector space V.
Df(r)(m)
The Fréchet derivative of an operator f at point applied to a perturbation
B SAMPLE COMPLEXITY IN CEMENT LEARNING
In this section, we provide an explanation sample complexity reinforcement learning. In our setting, the next-state information required by the reinforcement learning process is provided by PDE-based physical space. We assume that this physical space satisfies the property of probabilistic determinability of the future state, meaning that, as long repeated experiments sufficiently performed, the next state can be identified with probability at least one. We present a detailed explanation of the probabilistic determinability of the future state.
13

---

Published as a conference paper at ICLR 2026
B.1 DIRECTIONAL PERTURBATION NONLINEAR DYNAMICAL SYSTEM
To understand uncertainty propagation in complex dynamical systems, consider a nonlinear system governed by PDE-based evolution operator [yo, a), where yo the initial state and denotes system parameters or external inputs. We focus on a single dimension (or component) of the system state, say y(), and analyze how small perturbations in the initial condition affect the predicted future value vi (1) Suppose that we do not observe Yo exactly, but only perturbed version Yo + Ayo- where | I Ay reflects measurement uncertainty or modeling noise. natural question arises:
Whut the possible runge of future values Vi (t) under this uncertainty?
conventional linearized approximation assumes that the system's response to perturbations symmetric and smooth around yoo Under this assumption, we may apply first-order Fréchet
expansion
f + Ayo. a) ze f(yo. a) (4y0))A
and estimate the extreme values via projection onto the steepest ascent and descent directions. However, for highly nonlinear systems, this linear treatment is i when the system exhibits directional sensitivity, i.c.,
fwo 7) a) f'Vo a) + f(yo 7, a) f a).
To capture such asymmetry and nonlinear amplification, we adopt variational formulation. Let the system's response perturbation path 77 be described by a nonlinear variational functional:
Imnl / C(yo + sm, f). ss,
where L is a nonlinear Lagrangian that encodes the local effect of the perturbation path 7) and its derivative ii over time or space. The structural asymmetry of the system implies: I[n] & Z[-n]. highlighting that perturbations in opposite directions produce asymmetric responses.
We define the extremal range of the future value Vi (4)
the solution the following constrained
nonlinear variational optimization:
(Sm F()) F No. ((x), Vn((x) dp(x))
max
(10) + sup
(B.1)
mimin
(10 + inf B
{/1 (T()) You n(x), Vn(zx] du(z))
where F() denotes the nonlinear directional sensitivity of the system output with respect to the i-th component, and μ(z) is measure over the spatial domain M.
Remark B.1. While this formulation is presented for single component y()) it naturally generalizes to the full system state vector. This variational perspective offers a principled and scalable method for quantifying directional uncertainty in nonlinear, asymmetric dynamical systems, with broad applicability to PDE-based forecasting, control, and decision-making
As illustrated in Figure in the ODE case, the three blue points represent the lower bound, true value, and upper bound of the observed state respectively, resulting from observation noise. Since substituting different observed states into control equation yields different next-step predicted states, the three purple points correspond to the lower bound, true value, and upper bound of the next-step state.
We define prediction to be correct if it falls into the same grid cell the true next-step state. This correct cell represented by the purple box. Predictions that fall into other grid cells are considered incorrect.
The PDE example in Figure B.
illustrates how observation noise the initial state propagates to the next time step. Here, show the evolution over single time step. The red dashed rectangle represents the range of the observed initial state induced by observation noise. The use of a rectangular region stems from modeling the perturbation õy using the T norm, which captures worst-case
14

---

Published as a conference paper at ICLR 2026
ODE
PDE
(Ao)
y1 max
f((Ko. Vo) 40-m
1( a)
Nearest competing cell
Yo+
Correct cell
Y0
Yo by
Correct céjl
(x),)1)
min ff (Xo:Jo) + +yyaa)
Fo 1.
Figure B.1: A visual explanation of the probabilistic determinability of the future state.
uncertainty across all spatial dimensions. Under the assumption that the PDE solution operator is continuous with respect to the initial state in an appropriate Banach space, the open mapping theorem ensures that the reachable set of next-step states is also topologically open and smooth manifold segment.
To further characterize the image of this noisy initial set under PDE evolution, we adopt variational framework. Specifically, we consider a scalar output functional Fly] and examine extremal values under admissible perturbations dy U, where U C. M is a local uncertainty set defined on the tangent space of the solution manifold at yoo
Let the PDE solution operator be denoted v(t) S.(Mo). where yo is the initial condition and y() the state at time t. For a scalar output functional Fly) (e.g., wave peak or energy), our goal is to evaluate how F(S [Bol] varies under perturbations öy U C T M, M, where U is norn-bounded
uncertainty set on the tangent space of the state manifold at
Formally, the extremal values of the output functional under admissible perturbations given by the variational optimization:
Fmax = sup [B + 8y]]. 5yEM
Funi
F
inf SyEul
FlS + Sy]]
This general variational formulation allows for full nonlinear and directional sensitivity to be incor- porated. However, derive computable approximations and gain insight into asymmetry, we next truncate the variation up to second order.
As an illustrative example, consider the oneddimensional wave equation:
0,yy e² L.y.
y(z, 0) vo(x), 0,y(I, 0) 0.
Let the output functional be the solution at a fixed point and time: Fly] = y(0", At). We expand the perturbed output Flyo + dy] second order:
Flyo + 5y] = F[yo]= + // K(z)6y(o) 1 / 10.11=
H(x, af (8)x))))))) dx dx' + +(1by]))*
where K(x) is the first-order influence kernel, and Hxx is the second-orded interaction kernel.
The first-order term linear and sign-symmetric: it leads to the same magnitude of change whether by is positive negative. However, the second-order term breaks this symmetry. In particular, it captures nonlinear amplification and asymmetry: Flyo + öy] Flyo] * (Fly ōy] Flyol).
because the second-order term remains positive under sign reversal.
15

---

Published as a conference paper at ICLR 2026
Under the constraint [ S E, then compute bounds:
{ (() + // H(x, (xxxx))80 (o)(()))y(() do da' Fma = Flunl + sup oo Se
(B.2)
{ K(o)bolt) de lf H(z, do de' Se
{
Fuin Flyol +
inf
This example demonstrates how second-order variational expansion serves as concrete realization of the general nonlinear variational framework. It explicitly captures the directional asymmetry and curvature of the functional response to perturbations, which are essential features of PDE-governed systems with uncertain initial conditions.
B.2 SAMPLING SCHEME FOR IDENTIFYING THE TRUE DISCRETE CELII
In this section, we analyze how sampling strategy can be designed to identify the correct discrete cell that the system state should evolve into, and how the scale of such a cell should be chosen.
In real-world systems, perturbations Ay typically originate from measurement noise or numerical rounding errors, and are therefore inevitable. Nevertheless, it is reasonable to assume that the true state lies within a neighborhood of size Ay/2 around the observed value. Based this fact, we propose a local sampling strategy: perturbations are generated around the observed value within this neighborhood, and their evolution over one time step is computed via numerical PDE solvers. The resulting predicted outcomes are then statistically analyzed to infer the most probable cell in which the true state lies.
In the ideal noiseless case, such sampling scheme causes predicted outcomes to concentrate within the correct cell. As the number samples increases, the probability of correctly identifying the target cell converges to one.
This idea naturally supports the integration of physical knowledge into reinforcement learning. By repeatedly sampling perturbations around the initial state and propagating them through the physical model, we can estimate the state transition probability matrix with high confidence. This effectively transforms the reinforcement learning problem from model-free setting to a dynamic programming problem with known transitions. contrast, without such sampling, the system dynamics remain unknown, and exploration-basee algorithms must be employed to estimate both the model and the optimal policy.
To ensure probabilistic determinability of the future state, impose the following criterion: the number of predicted samples that fall into the correct cell must exceed that of any incorrect cell. Specifically, if the predicted uncertainty spans the correct cell and its two immediate neighbors, and if the correct cell is slightly larger in volume than its neighbors, then the correct cell can be statistically identified with high confidence. However, if one of the neighboring cells (e.g.. the upper cell) is significantly larger than the others, then it becomes necessary to ensure that the number of predicted samples falling into the correct cell still exceeds those falling into the largest competing cell. This places stricter constraints on the acceptable error magnitude Ay, or alternatively, requires more samples to suppress statistical variance.
To improve sampling efficiency, we further propose a variationally-derivee optimal importance sampling strategy. When the goal is to estimate the expected value of a functional output F + oy). importance sampling theory shows that the optimal sampling distribution p" (&y) should satisfy (ōy) & [F + 8y)1 p(5y), where p(öy) is the base distribution of uncertainty. In the small- perturbation regime, expand Flyo + 5y) using second-order Taylor approximation:
1 5y Hōy,
F(yo + 5y) F(0)) (VF, 8y) +
2
leading to the approximate optimal sampling density (ōy) x (VF 1,yy Hoyl Accord- ingly, we define sampling weight function: w(oy) := (VF, σy) + b H6y). which can be used to guide sampling toward directions with the most sigmificant influence on the predicted outcome. This minimizes estimation variance and achieves optimal sampling efficiency. The resulting scheme
16

---

Published as a conference paper at ICLR 2026
can be interpreted as physics-informed informed importance sampling strategy, grounded in variational sensitivity analysis and importance sampling theory.
We now quantify how many samples needed to reliably identify the correct discrete cell when using the optimal importance sampling strategy described previously. Specifically, we define the dominant cell as the cell that receives the largest number of predicted samples under given sampling scheme. Let C* denote the true (correct) cell and Cj denote any competing cell. Let nr" == (y +oy) E C*] and jyj [P[[(10 +dy) Cý}, where Py (ōy) Fyoo + +dy) p(ōy) the optimal importance sampling distribution.
Theorem B.1 (Sample Complexity for Correct Cell Dominance), Let K be the total number of cells, and let 8 (0, 1) be the desired confidence level. Assume that A == n²" Tmax > O, where Tmax = m Tz. Then, in order for the correct cell C to receive the largest number of samples with probability at least 1 it suffices to sample:
2(K -1)
N
A llg
8
samples from the optimal importance sampling distribution p" (ōy).
Proof. Let nc- ~ Binomial (N, x²) and n2 ~ Binomial (N, TTy) for each j & C To ensure that ne- > njj we use union bound and Chernoff-type large deviation inequality. The error probability that any incorrect cell receives more samples than C is bounded as:
N
P ( + C* nj 2 no-) S 2(K 1) exp
Setting the right-hand side less than equal to 8 yields the desired bound on N.
This result shows that the sample complexity scales inversely with the square of the gap A n² Tmax between the correct cell and the nearest competitor, and logarithmically with the number of cells and the inverse failure probability. This quantifies how distinguishable the correct cell is under the importance sampling distribution. The smaller the margin A, the more samples needed to overcome statistical uncertainty.
However, to preserve generality and narrative consistency, adopt uniform sampling strategy in the main analysis to estimate sample complexity, which differs from the optimal importance sampling scheme only by a constant-factor overhead.
B.3 PROBLEM FORMULATION
We consider two distinct scenarios concerning the availability of transition dynamics. These two cases correspond to two different learning regimes, each with its own complexity characteristics.
In the following analysis, we quantify the sample complexity under both settings. Based Lemma B. and Theorem B.2 we derive two different orders of sample complexity. The bound in
Lemma B. captures the general case without prior knowledge of the transition dynamics. In contrast, Theorem B establishes the bound under our setting, where the transition matrix is known through repeated sampling from the PDE-based physical space.
Therefore, if no sampling is performed in the physical space, the problem remains a standard MDP, and exploration is required learn the model. However, if sufficient sampling is carried out in the PDE domain, the problem effectively becomes dynamic programming task, and the sample complexity is driven by the cost of interaction with the physical model, resulting In the main body of the paper, we adopt the sample complexity bound given by Theorem B.2 In the following, we present several standard assumptions commonly used in the conventional MDP framework.
tighter bound.
Markov Decision Problems We consider episodic reinforcement learning in finite-horizon MDPs defined by < S, P. R. H >, where S and A are the finite sets of states and actions, P is the state transition distribution, the function R is a valuud function which is deterministic and belongs to the interval [0, 1). and the horizon H is the length of the episode. We denote by P(-|2r, a) and R(x, a)
17

---

Published as a conference paper at ICLR 2026
the probability distribution over the next state and the immediate reward of taking action state æ, respectively.
The agent interacts with the environment in sequence of episodes. Specifically, the agent follows a policy π S x [H] + A, which maps each state and time step to an action. The value function VI. S - IR represents the expected cumulative reward from step h to H under policy я, starting from state æ at time h. Under the previously defined MDP setting, there always exists an optimal policy x² attaining the optimal value function V (x) _ sup, V/," (3) for all x E S and h E [M].
Each policy π induces a transition kernel P7 (z/z) A P(ylæ, (æx h)) and reward function (x) A
For any function V S R, we define the transition operators:
(z) e Eyes Pf (vll)) The Bellman operator for policy is defined as (x)) & (=) + (z).
We assume that the reward function is known to the agent but the transition kernel is unknown (Under the no-sampling setting). The question we study how many episodes does learning agent follow policy ı that is not coppimml,, i.e., VM Vis, with probability at least 1 for any chosen accuracy c and failure probability 6.
B.4 FULL ANALYSIS OF SAMPLE COMPLEXITY BOUND
In this section, we systematically analyze the PAC sample complexity of reinforcement learning under two different assumptions regarding the transition probability matrix. We begin by examining the classical setting where the transition dynamics are entirely unknown. Lemma B.I provides baseline PAC bound in this setting, and reveal how the total sample complexity arises from key error components, including estimation error, exploration bonuses, and confidence adjustments.
Building on this foundation, we then analyze how our proposed setting-whig additional structure or prior knowledge is introduced modiiies the learning process. In particular, Lemmas B.2 and B.3 show how specific components in the original analysis are tightened or avoided, leading to an improved sample complexity bound.
Theorem B.2 (PAC Bound under our setting). Assume the transition dynamics of the MDP are known up to statistical confidence level 6, reducing the reinforcement learning problem to dynamic programming. Let the correct classification cell have predicted frequency p Ay)5), and let the
most probable competing cell have frequency q = Pmax- Then the number of forward predictions
(j)
required identify the correct classification cell with probability least 1 ô satisfies
log(1/5)
n ==
HSA
min, (Ayl)) - Pinax )
(5)
2
(SSA log (1)(1)))
Furthermore, we derive refined bound on the sample complexity n.
We also present in Lemma B.1 the PAC bound of the classical UCB-VI CCB-II&e algorithm under the setting where the state transition probability matrix is entirely unknown. The detailed algorithmic procedure of UCB-VI can be found in (Azar al. (017).
Lemma B.1. For any e, 8 < 1, the following holds. With probability at least 1 - ô, the algorithm
UCB VI produces a sequence of policies n² that yield at most O
SAAH 2a
log 1 episodes where
( VitH(so) (50) > e. This results in a total runtime for sampling of SAH log
Proof. Lemma B.I characterizes the PAC bound of the UCB-VI algorithm, result that has been rigorously established and thoroughly analyzed in prior work (see Azar et al. (2017)). To facilitate our later analysis, provide concise overview of the proof from Azar et al. (2017). with particular emphasis on the error decomposition framework and the core techniques used to bound individual error components.
Let V = (VE, 2 Vr. Vk. h denote the event where all computed value estimates are optimistic upper bounds of the optimal value function. Azar al. 2017) show that V holds with probability
18

---

Published as a conference paper at ICLR 2026
at least ò under standard induction arguments with appropriate bonus functions. Define the approximation emrors Wk,h Wa. (Zk,A) and Wk,h W where We. Ví V/ and Wa.A- Va h Then, under the event V, the following error decomposition holds:
1 Σ
(MRi 2VLMks + bonuss,; + FEs)) )
Wk.h S Wk.h =
WEAAI(ERAA-)) is martingale difference term. where MA.h- aax Wk. (E.A-
Ma is the quantity resulting from applying Bernstein's inequality to control the term (IA) and is also martingale difference, bomuss,; represents the bonus Vi++1) -
[PP PP (V+++
function for UCB-VI V1 at step under episode k. Ek def (P - p V (# is the estimation error of the optimal value function at the next state. addition, L here is a logarithmic factor depending on H, S, A, and the confidence level 8, which does not affect the polynomial complexity of the bound.
Based on the above statement, observe that the JCB-VI algorithm decomposes the total error into four components: My. MA.;. bonusg.i. and Eg.i- In Table B.2 we summarize the order of magnitude and the bounding techniques for each of these four error terms. The detailed proof process can be found in (Azar et al.) 2017). Therefore, the cumulative value error over all episodes and timesteps satisfies ΣR,h WK.h = O (HLYSAT) where collected across K episodes, each of length H. This implies:
KH is the total number steps
Σ Wk.h 0 = = = [HSI2LSSA)]
Suppose in the worst-case scenario, all of the first K episodes are not ج oplimal before an E - optimal policy learned, then we require:
H*L 62
SAH 62
K
S
Σ WR.h = =0 (1312 k.h
O SAKbad) = Kiad 50
0
log
To derive the total sample complexity, we note that sampling one episode and updating the respective variables has O(H) runtime. Thus, the total runtime across all Kbod episodes is O (H Kbad)
O SAH log }) This establishes the total sample complexity required to ensure that the learned policy is e-optimal with high probability.
Table B.2: Error decomposition in UCB-VI and bounding techniques
Error Term
Magnitude
Bounding Technique
Σx, Ma,
OHHTLL)
Azuma's inequality (martingale bound)
Ev,h Ma.
O(VTL)
Azuma's inequality (martingale bound)
En. Ex
OHHLSAT))
Chermoff-Hoeffding inequality pigeonhole principle
EA. bomusg,
O(HLVSAT)
Optimistic bonus design pigconhole principle
Remark B.2. Based Table B.2 have essentially clarified the analytical procedure through which the UCB-VI algorithu derives its sample complexity bound. Based on this, we highlight the key difference between setting and UCB-VI, and the core improvement introduced in approach.
UCB-V characterizes sample complexity by bounding the discrepancy between the empirical transi- tion probability matrix obtained from sampling and the true transition dynamics.
19

---

Published as a conference paper at ICLR 2026
In contrast, setting introduces notion of the probabilistic determinabilityy of state, which allows us to exploit multiple samples from the PDE-based physical space to ensure that the transition probability matrix in the RL formulation is fully known.
As a result, the increase sample size in our framework pertains to the sampling process in the PDE-based physical space. The complexity of this sampling process then be transferred the RL domain, ultimately yielding a refined bound on sample complexity.
In Lemma B.2 we provide the sampling requirement needed to resolve the correct cell with probability least 1 d.
Lemma B.2 (Sample Complexity for 5-Confidence Classification). Let the cell of correct classifica- tion have frequeccy p Ayd", and let the most probable incorrect (competing) cell have frequency q = Punx- To resolve the correct cell with confidence at least 6, the number of forward predictions required satisfies
log (L/() miny (ny')
O
(3) -maxx
Proof. Without loss of generality, we consider a randomly selected dimension j as an example.Le Ayli) be the probability mass of the correct cell in dimension j. and let Pmax be the maximum probability mass of any competing cell.
Suppose we sample n ii.d. trajectories, and let P denotes the empirical frequencies of the correct cell. We wish to ensure that with high probability (at least 1 - 6), pl lp-9l that the correct classification cell selected.
Let
Applying Hoeffding's inequality to gives:
Pr (10 2 () S 2exp(-2n22)
2nnc) )= = Pr 10
< 2exp(
2n( 1
2
2
Then enforce 2 exp
(-2n p
S
2log(2/6) [p-q]²
This gives the number of samples needed to distinguish between the correct and competing cells with probability at least To guarantee correct classification across all dimensions j. we must
take the worst-case (smallest gap) over all j. yielding: n min;
log(1/8)
(1) as stated.
Pinxx
Furthermore, Lemma B.3 provides the order of the total sample complexity under our setting.
Lemma B.3 (Refined Sample Complexity under Known Transitions). Assume the transition dynamics of the MDP are known up to statistical confidence level 8. reducing the reinforcement learning problem to dynamic progranming. Let the correct classification cell have predicted frequency p = Ay
and let the most probable competing cell have frequency q Pmnax Then the number of forward
(j)
predictions required identify the correct classification cell with probability at least 1 satisfies
log(1/8)
n ==
H,S.A, miny (Ay(0) - Pinax
2
-pillxx
Furthermore, by applying the p K relationship discussed in the main text, which states that Apli) ~ p 1 H. 5) we derive refined bound on the sample complexity:
O (HHS,A, log(1/6))
20

---

Published as a conference paper at ICLR 2026
Proof. Without loss of generality, consider randomly selected dimension j. Let p = Ay(0) be the predicted frequency of the correct classification cell, and let q = pidd denote the predicted frequency of the most likely competing cell in dimension j.
Since the transition matrix is assumed to be known up to confidence the reinforcement learning task reduces to dynamic programming problem. Assume we sample independent forward trajectories and let 9 be the empirical frequency of the correct cell. confidently resolve the correct classification, we want the empirical estimate to satisfy 1 17-pp 10-al with probability at least - 8.
Let 1p_ According to Lemma B.2 we have n 2 2 (p-g)²
At each step in the reinforcement learning process, every state-action pair requires the number of samples derived above guarantee correct classification with high confidence. Applying union bound over all such points yields a total sample complexity of:
=
0 (H.,,,,,
log(1/8) minj - r)2
Furthermore, have: p- 0 H. thus
O(H,). Substituting this into the (--
earlier bound yields: n = 0 (HRS,,, A, AoogAA.ooo///000 1og(1/8)) as claimed.
Remark B.3. Based on Lemma B. and Theorem B.2 we derive two different orders of sample complexity. However, we also leave open the possibility for further refinement. In particular, building upon the analytical framework presented above, future research may explore combining our treatment of transition uncertainty with recent developments on PAC bounds in robust MDPs, potentially yielding tighter complexity guarantees.
C DETAILED ANALYSIS OF THE FOUR REPRESENTATIVE A14S-RL ENVIRONMENTS
In this section, building on the theoretical framework proposed in the main text, we provide a detailed analysis of how surrogate discretization, spectral response, and reinforcement learning precision jointly affect sample complexity across four representative A14S-RL systems: (1) Tokamak plasma control, (2) Turbulent airfoil regulation, (3) Teppanyaki heat sequencing. and (4) Cart-pole stabilization. The first three environments are PDE-based systems, while the fourth represents typical ODE-based system. Notably, the first two environments are characterized as systems with strong boundary observability. Therefore, by conducting a detailed discussion of these four environments, we analyze how different environmental characteristics and underlying equation types influence the outcomes and properties of optimal resolution matching and computational cost, thereby enhancing the generality of our conclusions.
C.I THEORETICAL PRELIMINARIES AND PROOFS
Before proceeding to the detailed analysis of each environment, this subsection provides several general settings and supplementary explanations for key lemmas and conclusions referenced the main text. In addition, DoppendixDprroviiii collection of classical results from functional analysis and PDE theory that underpin the error analysis.
C.1.1 PROOF OF THEOREM
Theorem C.3 (Sample Complexity for 5-Confidence Classification). Consider repeated forward predictions from perturbed initial states, where the predicted frequency is the empirical probability of the next state falling into a given grid cell. Let p denote the predicted frequency with which the true next state falls into the correct grid cell, and let q = Pruxx denote the maximum predicted frequency among all competing cells. To resolve the correct cell with confidence at least 5, the number of
21

---

Published as a conference paper at ICLR 2026
forward predictions required is bounded by
1og(1/6))
n =
(C.1)
minj Ayy')
(i) Pmax
Proof. We proceed under the fundamental assumption of state distinguishability; which ensures that the propagation of imitial uncertainty is sufficiently bounded. This allows for the correct classification cell be identified, as repeated sampling from the input perturbation ball will generate predicted states that form distinct cluster primarily within that cell.
Consider independent forward predictions, each obtained by sampling a perturbation Ayo uniformly from the ball All I Ay} and computing ý1 f3o Ayo, a). Let X; be the indicator random
if Ø) falls in the correct cell
variable for the i-th prediction falling into the correct cell: X; 0 otherwise
The empirical frequency of the correct cell is: p + - X, By construction, E p and X; E [0, 1].
Since X2, .n,, are independent and bounded random variables with X; E [0, 1]. Hoeffding's inequality states: Pr(p pl 2 E) S exxp(-2ne" for any e > To distinguish the correct cell
from competing cells, we need ý to be closer to than to q. This is achieved if: 1 - pl < p-
Setting c 1p-gl p-a (since p q by assumption), we get:
2 exp -2n P PP- 2
p-q 9 Pr p p
2
2
For the classification to succeed with confidence at least
1- we require:
2
p 9
P 9
Pr pl
S 2 exp
-2n
5
2
Taking logarithms and solving for n: n 2 need to ensure correct classification across all dimensions j. Taking the worst-case scenario (smallest gap):
2
When multiple dimensions are involved, we
[p-q)=
2 2log(2/0)
2 log(2/ (7) min, (Ay()) -mmaa)
2 max (Ay()
G)
Pmmax)
In terms of asymptotic complexity (Big-O notation), the constant factors can be absorbed, yielding
log (1/(1/))
the final result: O miny (Ay()) p bound.
This completes the proof of the sample complexity
C.1.2 MODAL GROWTH RATE AND STATE SEPARABILITY
In AI4S-RL, agents interact with high-dimensional physical systems, where observations often contaminated with sensor noise or state perturbations. key concern whether such uncertainty grow through the system's nonlinear dynamics and cause classification ambiguity in the predicted future state.
To ensure robustness of the learning process, we must guarantee that small observation errors do not cause the predicted future state [ + Awo) to cross classification boundaries. This leads to the formulation of Lemma S.3 which characterizes sufficient conditions for achieving classification-ces cell separability..
Lemma S.3 (Fourier-Mode Separability via Fréchet Spectral Amplification). Let fe(y)) be the solution operator of nonlinear PDE with Fréchet derivative D)) acting Hilbert space H. Suppose the observation perturbation y) E H. with ||m]l < Ay. admits a Fourier expansion
n En He where (dx)) an orthonormal basis for H.
22

---

Published as a conference paper at ICLR 2026
Assume the Fréchet derivative propagates each mode (yv [ (t) (t) for some gain factor "/h (t) > 0. Then classification cell separability is guaranteed if: sup; (k(() < 1, which
implies (yo))il Ay. Consequently, the maximal admissible time step At satisfies:
At Cysstem inf 7/e (1).
Proof. Let = En MO E H be the perturbation with norm |,>/| I Ay, where {σr) } is an orthonormal basis. Since (fv( is linear, then we have:
Die(yo)lnl = Σ nm D. = E nu (t) vx(().
Now take the squared norm:
[Do) Σ n ma(1) (smppaact?)
Σ m sup |n11?
Therefore,
(De sup || s sup sup Ay.
&
k
To ensure (yo)i Ay, suffices that sup, (t) 1.
The spectral separability condition established in Lemma S.3 applies broadly to a wide class of PDE-governed systems. In general, when the system dynamics are governed by a nonlinear evolution equation, local behavior around state yo can be analyzed through Fréchet lincarization, yielding time-dependent linear operator D (z)()) By decomposing initial observation perturbations r) E H into orthonormal modal components and tracking the amplification of each mode under the action of Dfe(Uo), we obtain a spectral characterization of uncertainty propagation over time.
C.1.3 RELATING MODAL GAIN 7A (t) TO THE LEADING EIGENVALUE À1
In the framework discussed in the previous section, each mode experiences gain factor 7h (t), which reflects the sensitivity of the system to perturbations in that direction. When the maximal gain sup* 7 (t) remains below unity, the propagated perturbation stays within the classification boundary, ensuring separability and decision consistency. This condition naturally induces a wystem-dependert constraint temporal resolution: the RL decision step size must be sufficiently high to resolve the most unstable direction the spectral space. In practice, this analysis enables principled derivation of time-step bounds for variety of A14S-RL environments, thereby enabling robust policy learning grounded in physical dynamics.
Remark C.4. In nonlinear dynamical systems governed by PDEs, the rigorous characterization of forward perturbation growth should ideally be expressed in terms of the Fréchet-mode amplification factors (t), as introduced in Lemma S.3 These quantities measure the directional sensitivity of the solution operator Jt along each Fourier mode Φk. through the relation Dfe (Vooo))(o [áx] = mhk
where v/a (t) denotes the propagated mode shape. However, to simplify the analysis and maintain consistency across sections, we adopt the leading eigenvalue x₁ of the linearized system operator [(yo) as surrogate for modal growth, particularly in the computation of time resolution bounds and error propagation rates. The relationship between ds and mx(t) can be made explicit under the assumption of locally linearized evolution:
8yn (t) Aet OR
(t(() x edee
For sufficiently small t, Taylor expansion yields: (x(()) Rz 1 1+ det O(t?), so that sup, (t()) x )1. Therefore, replacing supk (t) with Aj in our main analysis introduces no change in asymptotic error scaling.
23

---

Published as a conference paper at ICLR 2026
C.1.4 PROOF THEOREM 2
Theorem C.4 (p-K Analysis for d-Dimensional Systems), For a d-dimensional PDE system, the forward projection enror rate is:
P p-1- 11/H, + 1 Ci/Ka
(C.2)
where the numerical error scales as Amoaal resolution limit:
CiK-KK " "y from spatial discretization. In the high-
p -= 0
(C.3)
H
Kd
Proof. The forward projection error rate, is defined as the fractional deviation from perfect prediction. Conceptually, this one minus the ratio the ideal state-cell size to the total spread of the predicted state under all sources of error. Let the ideal cell size be represented by the baseline observation uncertainty. Ay.
Ideal Cell Size
Ay
= =1 Total Error Spread
Total Error Spread
The total error spread is the linear superposition of three primary components identified in the text:
1. Intrinsic Error Growth (Aintrinsie): The amplification the initial uncertainty Ay by the system's dynamics over single RL time step At, H/H, This is governed by the leading eigenvalue Å1, yielding Amnrinsi À H1, Ay.
2. Baseline Observation Uncertainty (Aab):: The inherent noise floor of the system, which is Ay.
Numerical Surrogate Error (Amm)): The total error from discretization, Aomi, which theorem states scales as Atotal Ci Ky This scaling arises from the resolution matching conditions, where the dominant contribution to error that depends on the resolution ratio K is the spatial discretization of the underlying PDE in d-dimensional space.
Summing these components and canceling the Ay term yields the first result of the theorem:
H. +1+CK-
The high-resolutin limit is defined the case where H, and K + oo. To analyze the behavior of in this limit, we first rearrange the expression:
1
- H.
AI 1+
CCK-d
Let A + C,K-d In the high-resolution limit, H, - and K + oc, it is clear that - can therefore use the first-order Taylor series expansion for the function f(x) 1
H.
around r = 0. Substituting our expression for æ: 1 - pm 1 - (A म. +CKK- Then we obtain the approximate expression for p: p #.. K-d Since A and C1 constants for given system, the asymptotic behavior of captured by Big- notation
A
=
H
Kd
This completes the proof of the second part of the theorem.
24

---

Published as a conference paper at ICLR 2026
C.1.5 PROOF THEOREM 3
Theorem C.5 (Optimal Resolution with System-Dependent Scaling). For a physical system with slate space scaling S, K and action spuce scaling A, K under the computational balance condition H,S,A,, ~ Sw minimizing the computational cost
(108(1/0) E2
Cost(K)= H.KK
(C.4)
H. K
yields optimal resolution:
+ +++++2 2d
1/d
a 3 2d
1/a
K"
(ax (0[1-H111) H1:11
a+ 3
exp
dii
(C.5)
when 11, 2 11 >1
Proof. The total computational cost is given by:
Cost( K) H K log(1/0)
c2
H,
Kd
To find the minimum, we differentiate with respect to K. Since the prefactor H log(1/6) pendent of K, focus g()) K (1 H.
x² is inde-
1
To simplify the calculation,
Ka
logarithmic differentiation In (() (a + 3) In K + 21n
1
Hr
K
Differentiating both sides with respect to K: g'(K) 9(K]
atD
2d
At the optimal point
Katt
n. kd
K* we have g' ((K) = O, which gives:
a + B == 2d + +
(Ktya (1 IT,
1/d
Then (Ktd=
a++++2 (1+3)
Taking the d-th root K =
a+8+2d
(1+
We now analyze the expression for K in the limit where H, 2 A1 > 1. This implies that æ H a small positive value. We can rewrite the expression for K" as:
2d (2 x2222++++
1/d
K
(1-H,
a + 3
We focus on the second term, (1 -1- 1T Using the generalized binomial approximation, (1 - = xe 1 - ar for small æ, where x = H and a = -1/d: H1- se 1+ भे. For large \1, the argument the exponential small. Using the first-order Taylor approximation e 1 + y for small y: (1-H,") 1/d exp ) Substituting this back gives the final
approximate form for the optimal resolution:
2d aa++++ a B 1/d
3)
1
K"
exp
a
dy1
This completes the proof.
25

---

Published as a conference paper at ICLR 2026
C.2 TOKAMAK PLASMA CONTROL
C.2.1 ENVIRONMENT DESCRIPTION
In this task, we study magnetically confined plasma system involving the control of droplet-like plasma structure within a tokamak device. The objective is to manipulate the external magnetic field to maintain stable confinement of the plasma droplet and prevent wall contact (Degrave et al. 2022).
The governing dynamics are modeled by coupled magnetohydrodynamic system with deformable interface dynamics, comprising the following equations:
Op
+ V (pv) =0
(Continuity equation)
01
(
0v
+ V Vv
= - Vp + B
surface
(Momentum equation)
0t
ØB
V
(v B B
(Magnetic field control equation)
0t
06
+ V V = 0
(Deformation of the droplet's interface)
0t
OK
Ysurface
(Boundary conditions)
Here, p is the plasma density, v the velocity field, p the pressure, the magnetic field, and J = V x B is the current density. The term is the dynamic viscosity, is the magnetic resistivity, and F surfaaee represents interface forces. The function & level-set representation of the droplet interface, and its advection governs shape deformation. The final equation enforces curvature-related surface force balance, where a is the surface tension coefficient, is the local curvature, and Tsurface represents surface forcing.
The spatial domain discretized using a structured mesh aligned with the tokamak cross-section. A schematic of the cross-sectionss mesh layout is shown in Figure C.2 The state space includes grid-sam sam-le values of v, B, ø, and other derived quantities such as vorticity and magnetic energy. Control applied via discrete shape the confinement field. At each decision step, the reinforcement learning agent selects a control signal that reshapes the magnetic confinement field.
of actions adjusting time-varying coil currents Lumi (4), which
The reward function is defined as the negative of the minimum distance between the plasma droplet and the tokamak boundary, promoting stable suspension and avoiding wall contact. The long-term objective is to maintain the droplet within the desired confinement region while suppressing shape instabilities and maximizing equilibrium duration.
a. RL Space
b. Tokamak Device
PDE Space
Isoflux Line
Vessel
Plasma
Axis R,Z
position
Strike
points
Figure C.2: Structured mesh layout of the tokamak poloidal cross-section.
26

---

Published as a conference paper at ICLR 2026
C.2.2 ERROR COUPLING BETWEEN RL AND PDE SPACES IN A14S SYSTEMS
In this section, we analyze the sources and scales of error introduced in AI4S systems, especially the bidirectional projection between the RL space and the PDE-based physical space. As discussed in the main text and lemma S.3 for systems with observation noise, achieving an unbiased estimation of the RL transition kernel imposes lower bound on the temporal resolution of the RL system, which is determined by the intrinsic dynamics of the underlying physical process. Based on this property, we couple the errors in the RL space and the physical space, and ultimately derive an expression for the total system error. Then through p K analysis, we reveals the key factors that govern the overall numerical error the AI4S system.
We begin by introducing key notations to distinguish between two discretization structures:
The RL space is characterized by (H,,S, S,, A, denoting the RI. decision horizon, state discretization, and action granularity. The physical (PDE) space is represented by (H Sw Aw). corresponding to the physical evolution horizon, spatial resolution, and control parameter resolution in the underlying PDE model.
Without loss of generality, we assume that the spatial grid resolution satisfies Ar = Ay. Accordingly, we define Axp.um int the spatial grid spacing inside the PDE-based surrogate environment, and Axp.bd as the spatial spacing along the domain boundary. Similarly, Ax,, and Ax,,ba denote the interior and boundary grid resolutions within the RL environment, respectively.
Let the reinforcement learning policy operate on a coarse temporal scale, issuing actions every At, seconds, while the physical system evolves at much finer resolution Ate < At,. Consequently, each RL control interval consists of N, At,/ Ahw internal PDE integration steps. The total prediction error associated with single RL-AI4S step is then the cumulative result of local errors incurred at each fine-grained PDE time step.
We now analyze the composition of prediction error in AI4S systems. During each RL interaction cycle, the observed state must be projected from the RL space to the PDE space, evolved forward under the PDE dynamics, and then projected back into the RL space. This bi-directional mapping introduces two distinct projection errors, both of which may be amplified near domain boundaries due to the nonlinear characteristics of the governing equations.
Therefore, the total eтrоr consists of three main components: discretization error intrinsic to the RL space, numerica integration error accumulated within the PDE solver and coupling error caused by the nontrivial projections between RI. and PDE spaces.
Fine-Grained Surrogate Error in PDE-Based Surrogate Environments. We now analyze the surrogate error incurred at each fine-grained PDE step within the AI4S framework. This local error, denoted as A (R) arises from three primary sources: discretization in the spatial domain, imprecision in the control action, and interpolation near the physical boundary.
We consider the level-set PDE that governs the evolution of an interface function d: +v. Vo 0, where o(z, is the level-set function and v(z,t) is the velocity field. This PDE is discretized in time using a forward Fuler method and in space using central differences.
0e
of
onti-o" of
The numerical scheme used to evolve & at a single time step is: VndT the discrete approximation of the gradient at grid point x, and time tn, computed using central differences.
V, Vn = O, where
Alw
We now analyze the error introduced by this scheme, which is defined as the residual obtained by substituting the exact solution Φ(x, t) into the numerical scheme.
Temporal truncation error: Applying the Taylor expansion of p(x, t) in time at point i; gives:
do Ht
A12
0 n++ 1 = 6--h +
Ale
+ O(A)))
2
ae
The finite difference approximation of the time derivative becomes:
27

---

Published as a conference paper at ICLR 2026
0 +++1
0Φ m
Atw
86 dt² + (
+
Atw
0t
2
To establish the temporal truncation error, we consider the second-order Taylor expansion of the level-set function 622, t) in time at fixed spatial point z₁. In accordance with the physical model of tokamak plasma dynamics, the underlying equations exhibit parabolic-type behavior with smooth forcing and magnetic confinement fields. As such, we assume the solution Φ(I,, t) possesses sufficient temporal regularity.
In particular, assume & ( T) Q), which ensures that the second-order time derivative is uniformly bounded. That is, there exists a constant C' > 0 such that: * (z, 1) S C, V(x, () E Qr Under this regularity assumption, the finite difference approximation of the time derivative yields a local truncation error of the form: supremum of /
Time O(Atw), where the hidden constant depends on the
We consider the approximation of the spatial gradient Ve using central differences on a uniform grid with spacing Ar = Ay = h. For a function @ E H (52), we focus on the discrete gradient approximation defined as: 0 txe
Spatial truncation error:
d(#.-)) -1
2h
To evaluate the truncation error, we apply the first-order Taylor expansion o(#++1) = ((x,) + luy' (xx) - (12)). Substituting into the difference formula yields ((1441)=((((1))) ol d'(x)) O(h), which implies 0, de de
2h
+ O(h). A similar result holds for the y-component. Combining both 0z
components, obtain the discrete gradient V V + O(h), which holds pointwise under the assumed regularity.
In the context of PDE-based surrogate environments, the spatial truncation error depends the local grid region under consideration. Specifically, if we focus on the interior of the domain where Vo
the solution is sufficiently regular, the central difference approximation yields: V O Axp,,)) However, near the boundary or interface, the regularity of Φ is reduced due to physical or geometric discontinuities. In such cases, the interpolation error is governed by the trace theorem, and satisfies only: || 6b (a())) ((Ialxx.aall For further details the trace theorem and its implications for boundary error analysis, see Appendix D
Action-induced control error: In system, the velocity field v, which governs the evolution of the plasma interface, is not directly actuated but indirectly induced through a cascade of physical (MHD) equations. Specifically, the control input
interactions governed by the magnetohydrodyname are corresponds to the current applied to external boundary coils. The discretization of the action space introduces finite-resolution perturbations Adw, which in turn cause quantization-induce errors in the boundary actuation.
The control-to-interfface influence chain is follows:
Aaw
B 0n
J
V
x B -
The applied control signal perturbs the boundary magnetic field B Blan through electromagnetic models such as the Biot-Savart law or magnetic vector potentials 3B aa = O(Aaw). The magnetic field B evolves under the time-dependent induction equation dr v (v B) + V This is nonlinear parabolic PDE, where the diffusion term 17 2 dominates the short-time behavior. Over single control step, we can approximate the system as quasi-static and neglect the time derivative, resulting in - -1 V & V x (v x B). Under this elliptic approximation, boundary disturbances introduced by Au,, propagate into the interior via the smoothing effect of the Laplacian. By elliptic
regularity and the trace theorem, we obtain: 8B a = O(Aa
The momentum dynamics are governed by the MHD
equation:
+v. Vv = --Vp+JxB vv
surface:
0t
28

---

Published as a conference paper at ICLR 2026
where the Lorentz force J x B is the dominant actuator of plasma motion. Since J V B, both òB and 6J scale as ((Saw 2). This yields a force perturbation:
F Lorentz 6J x + x oB = O(Sa
and correspondingly sublinear response in velocity: Av O(Anu)))
Final per-step surrogate error: We now summarize the local surrogate error incurred in the Be + V = 0. Each term introduces a distinct
discretized interface evolution governed by: के approximation error: The time derivative θ,Φ' discretized by forward Euler, incurring O truncation enror. The spatial gradient V is computed via finite differences, contributing O O[]]Axpiaal in the interior and (I) (xp,xxl near boundaries. The velocity field v is perturbed by action discretization, with induced error O(Sa (22))
Substituting these into the PDE gives:
dn+1 6²
+ (v ((Aa
+01 (V)xO[]A
At
yielding the total local surrogate error per step
(k) = (Atw) + px. + O(Sa 1/2) + O(IAxx 113) + O(Aal/ ||Axp,
One-Step Prediction Error RL-Based Control Environments. In the AI4S control framework. the RL agent interacts with PDE-based simulator by issuing discretized action a, based observed state 8,, which itself a downsampled or filtered version of the true PDE state E V. The PDE system then evolves the dynamics over a time horizon At, and returns surrogate prediction error in this interaction arises from discretization of state and action inputs and integration-induced temporal accumulation.
state sy The
RL observation error as truncation error in PDE-Based environments: In PDE-based reinforce- ment learning environments, the agent does not observe the full continuous physical state y(x) e Rd but rather a discretized version 8,- E Rd obtained via sampling or interpolation over coarse spatial mesh. This introduces a spatial offset between the actual physical location z and the evaluation point ä used in downstream PDE computation: || E|| < Ax, |l. Such an offset mimics the numerical behavior of spatial truncation errors classical finite difference schemes.
Consider PDE used for control computation in the simulator: + v V = O, the derivative Ve evaluated not at the true point but at F. The following analysis parallels our earlier fine-grained surrogate error in PDE-based surrogate environments. from which we conclude that the truncation eтror induced by RL's observation of the state scales as Axp,ixx))
0e B8
Similarly, near the domain boundary an, the solution By the trace theorem, the restriction of ø to the boundary lies the fractional Sobolev space: ølan e (212). Consequently, interpolation errors on the boundary are governed by lower regularity, yielding I Φn (an) O(|Ax,,, bu 1/2), which mirrors the O(h rate of convergence found in boundary discretization analysis of finite difference schemes.
H' (92) typically lacks full smoothness.
In summary, the state observation errors in RL, when incorporated into the PDE solver, act implicit spatial discretization errors. Their scaling matches classical truncation theory:
[() Ax Il), 00( Ax,,,a 1/2). for boundary evaluation.
for interior evaluation,
! state
Action discretization error:
The agent selects action a,, from discretized control space, which governs the boundary actuation of the PDE system. The true control variable is approximated by a,, with finite resolution Aa,. This discretization introduces quantization emror in the applied boundary signal.
Similar to earlier fine-gra mme-grained surrogate error analysis in PDE-based environments, observe that this control mismatch propagates through the MHD system in a nonlinear and spatially diffused manner. The resulting perturbation in the velocity responilsee for evolving the system
29

---

Published as a conference paper at ICLR 2026
state-scales sublinearly with respect to the action resolution: Sa, dv ((aa'2). This velocity perturbation directly enters the PDE dynamics and contributes to the overall one-step prediction error the RL loop.
Temporal integration error:
During each RL step of duration At,, the PDE simulator executes N₁ / /t,/At,, At fine-grained steps with timestep Atw. At each substep, the forward Euler integration induces local truncation error:t; (m) = w()tw), so that the cumulative integration error over one RL
step becomes: Em Ti (n) O(A,,)
Summing all contributions, the one-step surrogate prediction error in RI. space satisfies:
ARL. C1 ,
C2ll C2 ball
+ CyAa
State (interior)
State (boundary)
Action resolation
Temporal propagation
Refined Total Error Decomposition with Time-Scale Separation. Based our previous analysis of surrogate error both the PDE and RL components, we now integrate all sources of discretization and numerical error into unified framework. Over a single RL step, the total prediction error can be decomposed as follows:
A Acotal
Cill Ax
+ C2 ||/2 +/2+ C3Da//2
+ CrAt
(interior space)
RL (boundary space)
RL (action space)
RI. (time)
Css
Air
(IAxp.imll + A/ + || 1/2
Dat/2 ||Axx, 11/2
(C.6)
Atw
PDE surrogute error (space action houndary), amplified by time scale separation
To ensure that the learned RL transition kernel remains robust under observation uncertainty Ay. we require that: A O(Ay). This constraint implies resolution matching condition across all discretization dimensions. Specifically, we obtain the following asymptotic scaling relations:
| ll ~ Axp,int || ~ Aa?/2 ~ Ao At, ~ Ay,
(C.7)
Aax
pbpp Ay² At Ay.
O Analysis. We now turn to the analysis of the error rate associated with predicting the next state AI4S environments based on an RL agent's s current observation and action. Specifically, we define prediction error event one where the next predicted state fails to fall within the correct discretized
RL state cell. Accordingly, we define the relative prediction error rate as: 1 Total Prediction where Ay is the size of the RL spatial grid. In a dinee-dimensional state space, it corresponds to the cube of the grid length.
Ay
The total prediction error in the system arises from a combination of three interrelated sources. First, observation noise introduces inherent uncertainty, stemming from either partial observability sensor imprecision within the reinforcement learning environment. Second, numerical surrogate error which quantified earlier in Equation equation C.6 the dissrreizztiioninnueee inaccuracies originating from both the RL and PDE components. Finally, the intrinsic growth of initial perturbations, as detailed in Lemma S.3 contributes additional error term due to the modal amplification behavior of the MHD system; specifically, this leads to where A denotes the leading eigenvalue and H, the RL plamning horizon.
uncertainty of order A Ay/ H,.,
Now we further quantify the numerical surrogate error based on the resolution matching relations previously derived. Given the grid refinement ratio K == Ax,/ AAxp which compares the spatial resolution of the PDE environment to that of the RL environment, and considering the combined error impact across three dimensions in state space, we obtain: Apoill & Ay) Substituting all three error components into the definition of p. we obtain:
Ay
p M Ay APOII
A1HH+ +1+
So in the limit of high-resolution settings, i.e., when H, and K are large, we obtain p =
O म. + Rz
30

---

Published as a conference paper at ICLR 2026
C.2.3 OPTIMAL COMPUTATIONAL COST ALLOCATION BETWEEN RL AND A14S
In this section, building on the previous analysis of error scaling, we provide detailed explanation of how to align the discretization parameters of reinforcement learning and A14S models to achieve optimal computational efficiency.
It is evident that higher fidelity in the physical model leads to more accurate predictions, which in turn accelerates convergence to the optimal policy RL. However, since AI4S simulators are typically executed under fixed computational budgets, it is crucial to ensure that the computational load in the RI. loop is of the same order as that of the AI4S module. Our objective to minimize total computational cost by optimally distributing resolution between the two systems.
According to Lemma B.I the sample complexity or computational cost of the RL component can be expressed COSTRL o HH
Meanwhile, the computational cost of the AI4S physical simulator is governed by its spatiotemporal resolution, and can be approximated by: CostA14S (.
To ensure computational consistency, we aim to align these quantities. This alignment is further guided by the theoretical result in Lemma S.3] where the RL horizon H, constrained by the intrinsic dynamics of the PDE system: H, Z 1 system.
X where x is the dominant growth rate in the linearized MHD
Based on the resolution matching conditions derived earlier, we assume the following scalings:
A, O(A.) 5, O(S.) (l H, H2/3 O(Tw)-
Here, Iw denotes the effective resolution of the A14S spatial grid. The scaling S, = (z arises from both the intrinsic dimensionality of the state space and the additional resolution requirements near the domain boundary. In tokamak plasma systems, the physical state space is inherently three- dimensional. Under uniform discretization, this implics that the number of internal state variables scales (zz)))
Besides, near the boundary a2, the PDE solution often exhibits reduced regularity-forman y H (Q)) then its restriction to the boundary lies in the fractional Sobolev space ( as stated by the trace theorem. This lower regularity implies that to maintain the same level of approximation accuracy at the boundary, the spatial grid must be refined further. In particular, the resolution must be squared in order to compensate for the smoothness loss, effectively contributing an additional factor of O (x))) from the boundary region.
The total RL cost then becomes:
COSTRI. H,S.A, 1og(1/0)
H
(C.8)
(Ay() Pmnax) prmax)2 )2
Hr
Our goal is to find the optimal refinement ratio K" that minimizes this cost expression. Formally, we solve the following optimization problem:
K
min KS0
C(K)
1
H.
K3
To find the stationary point, we differentiate C(K) with respect to K and set the derivative to zero:
4C dk
0. Denote
H. then C(K) K8
Taking logarithmic derivative:
KJ
d
2
Kd 3 Setting this derivative to zero and solving yields the optimal K*:
1/3
K*
So in the high-resolution regime where H, > 1, we approximate:
1(1- k;)
K 2
7
~
(Ze [
1//1
1/21
1/3
4(1
31

---

Published as a conference paper at ICLR 2026
C.3 TURBULENT AIRFOIL REGULATION
C.3.1 ENVIRONMENT DESCRIPTION
In this task, study unsteady aerodynamic system involving the control of a turbulent airfoil flow field. The objective to regulate the pitching angle of the airfoil to manipulate the surrounding flow structure and ultimately maximize lift (Portal-Porras et 2022).
The governing equations for the fluid dynamics are the two-dimensional incompressible Navier-Stokes equations:
du 01 + (u - - vAu
where u(z, y. (u(2, f), o( y. t)) is the velocity vector field.p(2,, y, denotes the pressure, is the kinematic viscosity, V is the gradient operator, and A is the Laplacian operator.
The system is also subject the incompressibility constraint: V еπоо sources and discretization structure that form the basis for our multi nuuti-resoluuiin analysis.
O. These equations define the
In this environment, we adopt a structured grid discretization for the spatial domain, employing an O-grid topology that conforms to the geometry of the airfoil. As illustrated in Figure C.3 a, the mesh progressively refined near the surface of the airfoil to accurately resolve the boundary layer. The state space is thus defined over this discretized grid. The action space is discretized into finite set of pitch angle increments, At each time step, the RL agent selects a discrete action corresponding to rotation of the airfoil about its chord line. The reward function is constructed based the instantancous lift coefficient, with the overall objective of maximizing its long-term average.
a. Airfoil
b. Teppanyaki
c. Cart-Pole
Interior Mesh
Airfoil
Angle
Boundary
Mesh
x
Position
Figure C.3: Visualization of Discretization Structures in Different Control Environments.
C.3.2 ERROR COUPLING BETWEEN RL AND PDE SPACES A14S SYSTEMS
Since the state space the turbulent airfoil flow field also features a division between interior and boundary regions, adopt analysis framework analogous that used for the tokamak device to analyze the composition of prediction error in A14S systems.
Fine-Grained Surrogate Error in Airfoil Flow Environments. We now analyze the surrogate error incurred at each fine-grained PDE step within A14S framework.In airfoil flow control tasks, the dynamics are governed by the two-dimensional incompressible Swiies-Stokes equations: d (u- V)u -- Vp + AAu. This PDE is discretized in time using forward Euler method and in space using central differences.
du
Temporal truncation error: Using the Taylor expansion time l,p write:
Ou
a
n+
A
+
O(A33
2
20
32

---

Published as a conference paper at ICLR 2026
Then the numerical derivative satisfies:
11 n+1 du m
u
Pu m dt +
Aew
At
öt
which implies a local truncation error of ((Atw) provided u E C?.
Spatial truncation error: For a function u E H (92), we focus on the numerical approximation of the convective term (u V)u and the diffusive term vÂu, which appear the Navier-Stokes momentum equation:
(u V)u Vp + vÂu. For the convective term, take one component
u(x, (), and use central differences to approximate: u
Bu s u(zi) u( -()==)) Since the "he
2h
function u (S2), the accuracy of the Taylor expansion is limited to O(1)). Applying Taylor
expansion, we obtain: u(2(41)-u(xx-) u (x;;) O(h), which leads to: u(x;) (((441)-a(ze)
2h
2h
+ (
For the viscous diffusion term, the Laplacian Au is approximated using second-order central and similarly in the y-direction. From Taylor expan-
differences: duu
42
sion: M( 1001:11)22uet ((l2 20(z, ux;))) + +h)). Hence, the Laplacian approximation also introduces an error of (h))
Therefore, both nonlinear convective and diffusive terms yield local discretization errors that can be bounded as:
(u V (u V)u + O(IIAxpimll))) A,u = Au + Axp.ix) I]).
On or near the airfoil surface an, the velocity field u may possess only H' regularity. By the trace theorem, its restriction belongs to: alan E (82) and the numerical interpolation error along the boundary grid satisfies: ||u - un (o2) = O[I 1/22).
Action-induced control error: In airfoil flow control, boundary actuation modifies wall velocities. Suppose this control is discretized with resolution Sog. The boundary perturbation in velocity is:
Sulag O( Because the viscous diffusion equation governs thort-time dynamics: du
vAu, we model the induced field internally by elliptic smoothing, yielding: Suln O(Aaw 1/2
Final surrogate error per step: Combining the above terms, the total local surrogate error per PDE step in the airfoil flow simulator is:
A (k) = O(Stw) + O(1 AX pmm + (/2) + (0Ilxxp
One-Step Prediction Error in RL-Controlled Airfoil Environments. Based on our earlier sur- rogate eтror decomposition for tokamak systems, we now analyze the one-step prediction error in RL-based control of unsteady airfoil flow. The underlying dynamics are governed by the 2D incompressible Navier-Stokes equations: ăt
du + (u V)u -p+ >Au. In this setting, the RL agent
selects an action a,- based on an observed state 8,, which a spatially discretized and possibly noisy version of the continuous flow field u. The PDE simulator uses a, to modify boundary conditions, evolves the fluid state over a macro timestep At,. and returns the next observed state s'.
RL observation error as fruncation error: Due to spatial discretization, the observation s, deviates from the true field u(a (z)) by at most the grid resolution: ||x - ã < ||,,,,, This offset behaves analogously to spatial truncation error in finite difference schemes. As previously analyzed, for smooth internal il int
state = O(I AX,, [). Near the airfoil surface, where boundary layers reduce regularity, trace theorem implies: ||u uz, ( = 0(|Ax,,,u|))
Action discretization error: The discretized control signal a, influences boundary actuation. Finite resolution Aa, introduces quantization error, which affects boundary velocity: óulan = O(Aa,). Via the viscous smoothing governed by vAu, this error diffuses into the interior: 8uls2 ((An//2)
Temporal integration error: Each RL step of length At, consists of N₂ = At,/Att fine-scale Fuler updates with local error O(Atw)) Total accumulation yields: ΣN, Tt (n) O(A)))
33

---

Published as a conference paper at ICLR 2026
Final surrogate error: Combining all contributions, obtain the total RL one-step prediction error:
C A1
+ C2 CarlAxx,hall
C3
CaAt
RI.
State (interior)
State (boundary)
Action resolution
Temporal integration
Refined Total Error Decomposition with Time-Scale Separation. Based on our prior analysis of surrogate errors in both PDE and RL components, now present unified decomposition of all numerical and discretization errors. Over a single RL step, the total prediction error can be expressed as:
A lotal
C 11 GxIAA,
C2114x,,,i 1//2+
Czdof/2
+ +Cadt,
RL (interior space)
RL. (boundary space)
(action space)
RT. (time)
+ Cs Air
1/2 + || p. 1/2 (1Axp.iml) + Aa
(C.9)
Atto
surrogate amplificd by time-seale separation
To maintain robustness of the RL transition kernel under observation uncertainty Ay. we impose the constraint: Awal O(ay). This leads to the following resolution matching conditions across
discretization components:
| || ~ AXp,int |l ~ An!/2 An ~ At, ~ Ay,
(C.10)
Ax,
Axp,bal Ay², At < Ay.
p-K Analysis. We now turn to the analysis of the error rate associated with predicting the next state AI4S environments based on an RL agent's current observation and action. We define the relative Ay
prediction error rate as: p = 1 Total Prediction Emerr where Ay is the size of the RL spatial grid. In two-dimensional state space, it corresponds to the square of the grid length. Similar to the previous analysis, the total prediction error comprises three key components: observation noise, numerical surrogate eπrоr, and intrinsic growth of initial perturbation.
Now we further quantify the numerical surrogate error based on the resolution matching relations previously derived. Considering the combined error impact across two dimensions in state space, we C. K Ay Substituting all three error components into the definition of we obtain:
obtain: Atotal
Ay
M + Ay + Apoaal p p=
A +1+ +
So in the limit of high-resolution settings, i.c., when H, and K large, obtain: p O + R2 The subsequent results on optimal computational cost allocation are summarized in Table of the main text. As this involves only straightforward calculations and scaling arguments, we omit detailed elaboration here.
C.4 TEPPANYAKI HEAT SEQUENCING
C.4.1 ENVIRONMENT DESCRIPTION
In this task, we study a two-dimensional heat transfer system representing a teppanyaki cooking surface. The goal is reinforcement learning control the heat distribution across the iron plate to ensure the fastest and most uniform cooking of food.
The thermal evolution of the plate is governed by the two-dimensional heat diffusion equation:
ÖT v.t) = Øt
PTT Oz
8ºT oy²
+ Q(x,
where T(x, y, t) denotes the temperature spatial position (z, y) and time t; is the thermal diffusivity coefficient, determined by the material of the iron plate; Q(x, y. represents the external heat source, which be actively controlled by the agent each grid cell (Koehler et al. 2024).
34

---

Published as a conference paper at ICLR 2026
The spatial domain is discretized using uniform Cartesian grid. Each grid cell corresponds to a controllable heating unit, allowing the RL agent to adjust the local heat input. As shown in Figure C.3 b. The action space is discretized into finite heating levels for each control unit. The overall objective is to optimize heating strategies that ensure efficient thermal convergence within the shortest possible time.
C.4.2 ERROR COUPLING BETWEEN RL AND PDE SPACES IN AI4S SYSTEMS
Unlike the previously introduced two environments, the teppanyaki setting adopts a uniform Cartesian grid discretization for the spatial domain. As result, there no distinction between interior and boundary regions, and no precision loss due to boundary regularity degradation. In the following, we adopt an analysis framework analogous to that used for the tokamak device to examine the composition of prediction enror in AI4S systems.
Fine-Grained Surrogate Error in Teppanyaki Heat Diffusion Environment. We now analyze the local surrogate error incurred in each fine- the-grained PDE step under the teppanyaki heat diffusion setting. The governing equation for the temperature field (xx y, t) over the iron plate is the two- dimensional heat diffusion equation: 2rt = a&T + (ax, y: I), where a the thermal diffusivity and Q is the external heating source, modulated by the control input.
Temporal truncation error: We discretize the temporal derivative using the forward Euler scheme:
T. Th = TA + QT, where T, is the discrete Laplacian at grid point i, Applying Ar.
Taylor expansion, the numerical derivative satisfies:
T Ty
ØT
OT
R + O(A)),
A
At
at
OF
k.j
i.j
which implies a local truncation error of ((A,,) provided T E C²2 We obtain the temporal truncation error:Tume O(Atw).
Spatial truncation error: Tit1, 2T. + TT-11 Ti.+1 - 2T; j + Tij-1
The Laplacian operator is approximated by central differences:
AT ~3
(Ar)²
(Ay)²
Assuming T (S2), the second-orded spatial derivatives exist in the weak sense, but due to This
limited regularity we expect that: A,T AT + O(IAAx)1)) This yields a spatial truncation error:
Tspace ([1 Axpll).
Action-induced control error: this environment, control is implemented by adjusting the external heat source term Q(x, y, 1) in the heat diffusion equation. Let the continuous control signal be denoted by a(I), and suppose it discretized with resolution Au,, into a piecewise constant control aw (t). Assuming that the mapping from the control input to the heat source is linear, we obtain: Q(x, y. = (( y. t) + 8Q(x, (z, y. t), where Q is the nominal source corresponding to (t), and the perturbation in the source is: 8Q(z,
(z. y, fōa) O(Daw)
Since the heat equation is linear in Q, the perturbation Q propagates linearly to the temperature field. Hence, the resulting temperature deviation induced by control discretization satisfies: 8T O (Aa,))
Total per-step surrogate error. Combining all error sources, the total surrogate error in the tempera- ture prediction per PDE step is: A ((Atw) OIIAx,,II) + O(Aaw).
One-Step Prediction Error in RL-Controlled Heat Diffuusion Environments. the teppanyaki plate task modeled by the heat diffusion equation, the RL. agent interacts with PDE simulator by issuing a discretized action a, and observing discretized temperature field 8,, which derived from the true temperature field Tæ y, t) through uniform spatial sampling. The simulator then evolves the system over a time step At, and returns an updated observation s/.. The one-step prediction error in this process arises from three main sources: state observation error, action discretization error, and temporal integration error.
RL observation error: The observed state s, is obtained by projecting the continuous temperature field T(x, y. t) onto a coarsely discretized uniform Cartesian grid. This introduces a spatial obser- vation error due to the offset between the true location (Ir, y) and its discrete representation ã,
35

---

Published as a conference paper at ICLR 2026
Assuming T E (Q), this interpolation sampling procedure introduces first-order truncation error: [2, - r) = A(,)Ax,,.mll)) I]). This matches classical interpolation theory over uniform grids and reflects the interior spatial error due to limited resolution in RL's perception of the state.
Action discretization error: The action space discretized with resolution Aa,. The true continuous control is approximated by its discretized counterpart o,,, and thus the perturbation in the heat source term (ax y. t) is: 5Q(z, y, t) = Q(a) ((a,) = O(Aa,). ). Since the heat diffusion equation is linear Q, this error translates directly a temperature deviation of the same order: 87 = (DSa,)
Temporal integration error: The PDE simulator advances the solution over a time horizon AI, using N₂ A Alw steps of the forward Euler scheme. Each substep incurs local truncation error: Tt (w yielding a total accumulated temporal error: >m=1 T O())))
Combining the above components, the total one-step prediction enror induced by RL interaction with the PDE simulator in the heat diffusion system is given by:
A RL C CillAx,,im
+
C2Sa, Action discretization
Cs
State observation
Temporal integration
This decomposition forms the basis for analyzing the interaction fidelity and prediction uncertainty of RL agents in physical systems governed by parabolic PDEs like heat diffusion.
Refined Total Error Decomposition with Time-Scale Separation. Based on our prior analysis of surrogate errors in both PDE and RL components, now present unified decomposition of all numerical and discretization errors. Over a single RL step, the total prediction error can be expressed as:
C₁ At, Atg (1Axpl + Daw)
Atotal C |x
C2Sa, (action discretization)
C3At,
RL (state space)
(temporal propagation)
PDE surrogate error, scaled by time resolution
(C.11)
To ensure the error stays below the acceptable uncertainty level Ay, require: Apotal O(Ay).
This yields the following matching constraints for all discretization parameters:
Axxll ~ || ~ Sa, Dow At,, ~ Ay, with Alw & Ay.
(C.
p-K Analysis. We now analyze the prediction error rate in the heat diffusion AI4S environment. We define a prediction eтror event as one in which the AI4S simulator's 's output state does not fall within the correct discretized cell of the RL state space. This gives rise relative prediction error Ay
Total Prediction Emor where Ay is the size of the RL spatial grid. In two-dimensional spatial domain (as in heat diffusion over a cooking surface), Ay corresponds to the square of the grid length. The total prediction error contains the following components: observation noise, numerical surrogate error and intrinsic error propagation. Over an RL planning horizon H,, the error grows diffusively at rate A₃. yielding an error amplification:2, Ay/ H,.
rate: p
Substituting all components into the definition of we obtain:
Ay
P =
A/ - By CAAyKK2 =1 A/HH +1 C1/K2'
So in the limit of high-resolution settings, i.e., when 11, and K are large, we obtain: p O O रス The subsequent results optimal computational cost allocation summarized in Table of the main text. As this involves only straightforward calculations and scaling arguments, omit detailed elaboration here.
36

---

Published as a conference paper at ICLR 2026
C.5 CART-POLE STABILIZATION
C.5.1 ENVIRONMENT DESCRIPTION
In this task, consider the classic inverted pendulum control problem, where the objective is to apply a horizontal force to a cart in order to maintain a pendulum in an upright and stable position.
The state of the system is represented by a 4-dimensional vector: x [x denotes the position of the cart, I the velocity of the cart, & the angle of the pendulum from the vertical (positive clockwise), and the angular velocity of the pendulum. The schematic of the discretization of position and angle the inverted pendulum environment is shown Figure C.3
0) where
The control input is a horizontal force u(t) E [-Fmax: F F.ax], which we discretize for reinforcement learning purposes. The agent selects actions from this discretized control space to influence the system dynamics. The physical parameters governing the system include: M: the mass of the cart, m: the mass of the pendulum, I: the length from the pivot to the center of mass of the pendulum, g: gravitational acceleration.
Using the Lagrangian formalism, the nonlinear dynamics of the system can be derived as (Nagendra)
al. 2017): :
9 sin 0 cos 0 [u+ siu +m
cos A M+m
u + md (2 sin B cos
3
M
These equations define the state transition model and control response of the inverted pendulum system. The reinforcement learning agent interacts with this environment by issuing discrete actions at fixed time intervals, aiming to stabilize the pendulum around the upright position. The reward function is typically defined penalize deviations from the vertical orientation.
C.5.2 ERROR COUPLING BETWEEN RL AND ODE SPACES IN AI4S SYSTEMS
Unlike the previous three environments, the control dynamics in the inverted pendulum environment are governed by an ordinary differential equation (ODE) system. But the error analysis procedure remains analogous to that used in the earlier cases. In the following, we analyze the composition of prediction error for the inverted pendulum environment.
Fine-Grained Surrogate Error in Inverted Pendulum Environments. We now analyze the surrogate eгror incurred each integration step in the inverted pendulum system. Since the dynamics are governed by ordinary differential equations (ODEs), the error sources are reduced to: temporal discretization, state observation through interpolation, and action quantization. We focus on the position variable xx() of the cart and quantify its one-step prediction eтror.
Temporal truncation error: The ODE describing the cart's horizontal motion is: ă f(x, u), where u E [-Fmax: Fmax] is the external force and x = [z. à 0. is the state vector.
Using forward Euler integration for both velocity and position: 4n+1 in + Alp x + i the local truncation error from Taylor expansion becomes: x n1+ x ==+1_=" Ar. 2n + 2
Aly
A12 3n O(At.)). Then the numerical derivative satisfies:
At O(At))
From the expression of F in the inverted pendulum dynamics, we know that Tn is bounded. Hence, we obtain the temporal truncation error: Time (At,)).
State diseretization error The state is often represented in coarsely discretized state space. Suppose the cart's 's true position & E IR is approximated by a discrete cell center ă. Using linear interpolation neireot projection, the induced emror is: Iz ãl OxI]xx,))) which propagates into the dynamics via the nonlinear function fxx u). Assuming Lipschitz continuity of f. this implies: A state O(IAAx
37

---

Published as a conference paper at ICLR 2026
Action-induced control error: The control force is discretized into finite set, and the applied value Uw approximates the true value with resolution Aa. The dynamics respond linearly to changes in u, from the equation for I, so: Acction O()a.).
Final per-step surrogate error: Summing all components, the total one-step surrogate error in predicting the cart position is: A (k) ((A),) + O(I)x, + O(Aaw).
One-Step Prediction Error in RL-BBase Inverted Pendulum Control. this task, we analyze the surrogate prediction error that arises when an RL agent interacts with the dynamical system of an inverted pendulum. The dynamics are governed by a nonlinear second-order ordinary differential equation (ODE), and focus our analysis on the position variable =:(), which affected by both control input and physical evolution. The prediction error consists of three major sources: temporal discretization, spatial resolution (observation) error, and control discretization error.
The ODE simulator advances the solution over a time horizon At, using N₂ At,/ At steps of the forward Euler scheme. Each substep incurs a local truncation eπror: T ()) = ( yielding total accumulated temporal error: Σ T (a) OAt,))
Temporal integration error:
Suppose the RL agent does not access the true continuous state attn), but instead receives a quantized observation (t,,). derived from uniform discretization over spatial grid resolution Ax,, In one-dimensional observation with linear interpolation or rounding. the state projection satisfies: (zr((n.. - =((.)) < Az,. This discrepancy is formally equivalent to spatial interpolation error and yields spatial state uncertainty: Tuate O(rr))
Observation-induuee spatial error:
Action discretization error: The control input u is issued from a discretized action set E Ar, such that: utx,,) = ū(tx)) + du(t,), du O ((aa)), where ū(t,) is the true continuous control signal and du is the discretization mismatch. From the expression of म. this control deviation contributes linearly to the evolution of position through: 8x = 5t,5 = Ar, áu
M++ = O(Aa,) We
treat time as a fixed constant when integrating the accumulated effect the action, rather than considering it as an error term.
Total error: Combining the above contributions from temporal, spatial, and control discretization errors, the one-step prediction error in the inverted pendulum environment satisfies:
! pendulum
CLAT State resollution
C2Dar Action resolution
C2At
Temporal intogration
Refined Total Error Decomposition with Time-Scale Separation. Based our previous analysis of surrogate error in both the ODE and RL components, we now integrate all sources of discretization and numerical error into unified framework. Over a single RL step, the total prediction error can be decomposed as follows:
At, Atg (1Axpl
Apotal
C A1IAx,11
C2Sa,
+
CsAt,
C₁
Daw)
RL (state space)
(action discretization)
(temporal propagation)
PDE surrogate error, scaled by time resolution
(C.13)
To ensure that the learned RL transition kernel remains robust under observation uncertainty Ay. we require that: Apoul O(Ay). This constraint implies a resolution matching condition across all discretization dimensions. Specifically, we obtain the following asymptotic scaling relations:
Ax,||
|
Aa,
Ddw At,
Ay, with Alw Ay.
(C.14)
p-K Analysis. We now analyze the prediction error rate in the cart-pole environment. The total prediction error contains the following components: observation noise, numerical surrogate error and intrinsic error propagation. Over an RI. planning horizon H,, the enror grows diffusively at rate X yielding an error amplification: X, Ay H,. Substituting all components into the definition of p. we obtain:
Ay
1
0=1 AAA + Au C,Ay/K = X/HI, +1+ C₁/K
38

---

Published as a conference paper at ICLR 2026
Since there is notion of state grids in the inverted pendulum environment and all state discretizations are first-order, the order of is one. So in the limit of high-resolution settings, i.c. when H, and
The subsequent results on optimal computational cost allocation are summarized in Table of the main text. As this involves only straightforward calculations and scaling arguments, we omit detailed elaboration here.
K are large, we obtain: p 0 7 +
C.5.3 EMPIRICAL COMPARISON OF TABULAR RL. AND DEEP Q-LEARNING
We designed an experiment to test our theoretical findings on different learning algorithms. Using data-driven inverted pendulum (CartPole), we compared the sample complexity of Tabular RL and Q-Learning. This comparison was performed under various discrete observation resolutions.
Experimental Setup The CartPole system consists of cart moving along a horizontal track with hinged pole on top. The agent's goal is to apply horizontal forces to maintain the pole in an upright position. We adopt a two-stage workflow: (i) seate-action-ext--see tuples are collected from physics- bassec simulator and used to train an MI.P-based dynamics model, and (ii) this learned model serves as the environment for RL training.
Reward Function To provide smooth learning signals and mitigate oscillations, employ a shaped reward function:
(1 a a000a 1 -|0.. if [0cl S 12" and |arel S 2.4
F1
(C.15)
otherwise
where Og and O, are the pole angle and angular velocity, respectively. The coefficients o, >> 0 penalize deviation from the vertical and excessive oscillation. An episode terminates with zero reward the pole falls or the cart leaves the track.
Learning Methods We train two agents under identical discretized observation resolutions:
Tabular RL: Implemented using the value iteration algorithm, with environment dynamics derived from physics- based simulation and interpolated the RL observation space, without neural network approximation.
Q-Learning: Implemented using Deep Q-Network (DQN) with discretized state and action spaces.
Performance measured by the number of samples required to achieve value function error below % and the corresponding wall-clock runtime single NVIDIA RTX 4090 GPU.
Results The detailed experimental results are presented in Table CS and Table C.4 The results show that Q-Leaming is significantly more sample-efficien than tabular RL, which is expected due to the generalization capability of the neural network. However, both methods show a clear sensitivity to the resolution parameters H, and K. Optimal performance is not achieved at the highest or lowest resolutions but rather at an intermediate point, which is consistent with our thcoretical analysis. For tabular RL, the optimal configuration is found at K 1.5 and log H, 1/3. For Q-Learning, the optimum shifts to K = 2.0 and log H, = 1/2, suggesting that DQN can better leverage a higher-fidelhe physical model. These results strongly support central claim that even for deep RL methods, the resolutions of the physical simulation and the RL algorithm must be carefully coordinated achieve optimal performance.
D FUNCTIONAL ANALYSIS TOOLS FOR PDE ERROR AND LIMIT ANALYSIS
In this section, we provide a collection of classical results from functional analysis and PDE theory that underpin the eтror analysis, convergence behavior, and design of discretization schemes in AI4S- RL environments governed by PDEs (Eaan 2022 2011). These tools help quantify how physical solution properties interact with discretized approxımations, especially under limited regularity, complex boundaries, and dynamic state evolution.
39

---

Published as a conference paper at ICLR 2026
Table C.3: Cart-Pole System: Sample Complexity (o))o(N) for Tabular RI. Achieve Value
Function Error 0.01.
logn H
K=1.0
K=1.5
K=2.0
K=2.5
1/8
4.51(4.58/4.45
4.50(4.53/4.46
4.51(4.57/4.44)
4.54(4.61/4.48)
1/6
4.26(4.32/4.19)
4.25(4.31/4.18)
4.26(4.30/4.23)
4.29(4.3//4.23)
1/5
4.10(4.17/4.04)
4.09(4.13/4.06)
4.10(4.16/4.03)
4.14(4.21/4.08)
1/4
3.91(3.98/3.85
3.90(3.94/3.87)
3.91(3.95/3.88
3.94(4.01/3.88)
1/3
3.77(3.81/3.74)
3.65(3.71/3.58)
3.74(380/33677
3.88(3.95/3.82
1/2
3,95(4.02/3.89)
3.83(3.87/3.80)
3.92(3.99/3866
4.06(4.13/3.99)
4.25(2.2
4.13(4.20/4.07)
4.22(4.26/4.19)
4.36(4.42/4.29)
2
4.55(4.62/4.49)
4.43(4.47/4.40)
4.52(4.59/4.46)
4.66(4.72/4.59)
Table C.4: Cart-Pole System: Sample Complexity (N) for Q-Learning to Achieve Value
Function Error < 0.01.
log H,
K=1.0
K=1.5
K=2.0
K=2.5
1/8
3.20(3.26/3.13)
2.95(3.01/2.88)
2.82(2.86/2.79)
2.90(2.97/2.844
1/6
3.05(3.11/2.98)
2.80(2.84/2.77)
2.72(2.78/2.655
2.80(2.86/2.73)
1/5
2.90(2.94/2.87)
2.75(2.81/2.68)
2.68(2.72/2.655)
2.75(2.82/2.99)
1/4
2.80(2.87/2.74)
2.70(2.74/2.67)
2.7((2.3/226611
2.72(2.76/2.69)
1/3
2.78(2.85/2.72)
2.69(2.73/2.666
2.66(2.72/2.59)
2.71(2.78/2.655)
1/2
2.65(2.69/2.22) 2.67(2.33/2.600) 2.75(2.79/2.72)
2.67(2.44/2.61)
I
2.85(2.91/2.78)
2.78(2.82/2.75)
2.78(2.85/2.72))
2.82(2.86/2.79)
2
2.95(2.9//2.22) 3.10(3.17/3.44)
2.93(2.88/2.77)
3.00(3.07/2.94)
D.1 MOTIVATION: WHY DISCRETE FORMATS AND FUNCTION SPACES MATTER IN AI4S-RL
In PDE-governed A14S-RL tasks, control and prediction rely the accurate approximation of continuous state evolution. However, neural networks or other function approximators inherently operate over discrete representations. Therefore, the erro between the true PDE solution and the learned approximation depends not only on the architecture but also on how the discretization aligns with the analytical properties of the solution. Tools from functional analysis allow to assess:
How regular or irregular the true solution is (c.g., near boundaries under shocks);
How the discretization grid or function basis must adapt to preserve convergence;
How stability and generalization are affected by boundary singularities and approximation limitations.
The following theorems provide essential insights into these questions and can be naturally organized along a logical chain used in theoretical PDE analysis:
Start with basic control over functions via gradients (Poincaré Inequality),
Link smoothness to representability and generalization (Sobolev Embedding),
Ensure PDE well-posedness (Lax-Milgram),
Connect numerical approximations to best approximants (Céa's Lemma),
Use compactness to show convergence under refinement Beeliih--Konrrchhyy
Extend to temporal-spatial regularity (Aubin-Lions),
Analyze boundary effects (Maximum Principle, Trace Theorem),
Control physical evolution through energy bounds (Energy Estimates).
The materials presented in this section well-established results in computational mathematics and are not original contributions by the authors. Our aim is to provide readers paaieers-pptticuarrlyy those interested in the error analysis of PDE-based surrogates-with introductory overview of the
40

---

Published as a conference paper at ICLR 2026
foundational knowledge in this field. We emphasize concepts essential for understanding boundary error behavior in PDEs, especially those involving the trace theorem, which is repeatedly applied throughout this paper.
D.2 CORE THEOREMS AND THEIR ROLES
1. Poincaré Inequality. Let H (S2) bounded domain & I Then:
| & Vl|
Interpretation:
This inequality ensures that if the gradient of function is small, then the function itself must also be small in the L sense. It captures how global behavior can be controlled through local derivatives. Use: In AI4S-RL systems, this inequality provides the foundation for energy- -ased bounds, helping ensure that learned control policies do not cause unbounded or physically unrealistic growth system states during training or rollout (Evans) 2022).
2. Sobolev Embedding Theorem. For sufficiently smooth domains and appropriate ( - LP()) or C* ($2).
Interpretation: This theorem allows one to infer integrability continuity from Sobolev space membership, linking smoothness to spatial regularity. Use: In AI4S-RL, this theorem helps justify the of neural approximators, showing under what conditions function learned in weak (Sobolev)
sense can be expected to generalize as a strong or continuous function (Evans 2022).
Lax-Milgram Theorem. Given bounded, coercive bilinear form a(;-) on Hilbert space V, the variational problem:
a(u, v) = f()), V
has a unique solution E V.
Interpretation: This theorem guarantees that weak formulations of PDEs are well-posed under mild assumptions. Use: It is critical for framing A14S-RI problems involving constrained learning, such solving inverse control problems governed by PDEs. Ensures the existence and uniqueness of target states (Evans) 2022).
4. Céa's Lemma. Let solve the variational problem and un V/, V be the approximation. Then:
C inf llu-valle V.EV,
1u-wullv- S
This lemma asserts that the best possible error of a numerical approximation is bounded proportionally by the best projection error. Use: In AI4S-RL, this guides the choice of rep- resentation spaces (e.g., basis functions, neural network architectures) and ensures that approximate policies value functions can converge quasi-optimally (Hou 2003).
Interpretation:
5. Keldach-Kondrachov Compactness Theorem. If 92 is bounded and u, is bounded in H+ then subsequence converges strongly in (2)))..
Interpretation: This compactness result ensures that bounded sequences of approximate solutions have strongly convergent subsequences. Use: It provides a theoretical foundation for showing convergence of learned policies or solutions in AI4S-RL as the discretization or model resolution is refined (Evans) 2022).
6. Aubin-Lions Lemma. Let No + X X1 be Banach spaces with compact embedding. If bounded in (2 Xo) and 0,u, is bounded in L* X1), then u, is relatively compact in (P T; X).
Interpretation: This lemma bridges temporal and spatial regularity, ensuring strong convergence when time derivatives and spatial norms controlled. Use: Vital in PDE-driven RL settings with temporal evolution. For instance, helps demonstrate that learned state sequences from policy iteration converge true continuous dynamics [Lions]]] 1969..
41

---

Published as a conference paper at ICLR 2026
7. Maximum Principle. Let solve elliptic parabolic PDE with suitable structure. Then:
max u
max t. 052
This principle ensures that solution extremums occur on boundaries, preserving physical constraints. Use: In AI4S-RL, this is crucial when dealing with bounded domains and enforcing reward shaping or physical limits. Helps constrain exploration near boundary zones (Majda) et al. 2002).
Interpretation:
Energy Estimates.
Multiply a PDE by u or Ozu, integrate over space-time, and derive:
E(t) = S/ (10,ul= Vul²) dr is bounded decaying.
These estimates show that system energy either remains bounded or decreases over time. Use: Critical for stability analysis of RI. agents physical systems. Ensures that learned
Interpretation:
dynamics don't 't lead to energy blow-up or unrealistic oscillations (E) 2011 Majda et al. (0022).
9. Trace Theorem. For u E H (S2), its restriction to the boundary 0Q satisfies:
(u)) := ulan E H1 (
Interpretation: Functions in Sobolev spaces lose half a derivative when restricted to boundaries. Use: This loss of smoothness explains why control near boundaries harder in AI4S-RI. systems and why special handling is needed in boundary-constrale reward or control definitions (Evans 2022).
D.3 LINKING THEORY TO PRACTICE
These tools collectively help answer: How should discretize state, design control inputs, and construct reward functionals when the solution is only partially smooth, or when boundary behavior dominates the error? Trace and Poincaré help identify when boundary refinement is necessary; Sobolev embedding and Rellich-Kondrachov justify the use of compact approximators; Lax-Milgram and Céa give the mathematical foundation for why value function approximators can converge. Understanding and incorporating these results into the AI4S-RL model design allows for principled control over generalization, robustness, and physical consistency in the learning loop.
E EXPERIMENTAL DETAILS: TEPPANYAKI THERMODYNAMIC CONTROL MODEL.
E.1 PROBLEM FORMULATION AND PHYSICAL MODEL
E.I.1 PHYSICAL SYSTEM DESCRIPTION
We construct a teppanyaki temperature control simulation based on the two-dimensional unsteady heat conduction partial differential equation. The computational domain a uniform square metal plate Ω [0, 1 x [0, 1] m². The system is equipped with H, = 2 power-adjustable heat sources. These sources are fixed at spatial positions x! (1) (0.25, 0.5) and x (0.75, 0.5). The control variables are the heating powers (t) E [0, 1000] W for each source.
Three dishes are placed on the plate. Each dish i modeled as a square region D, C Z with side
length La = 0.1 m. Their geometric centers are located at x = (0.4, 0.5), x = (0.5,
and x (3) (0.6, 0.5). The target temperatures are {70, 80, 90) "C, and the cooking durations are (50, 55, 60) seconds. The control objective is to minimize the average temperature deviation of all
dishes within their cooking windows over a total duration of Trotal 60 seconds.
The regional average temperature of dish at time t is defined as:
(t) = [D. /f.
1
T(x, t) dx
(E.I)
42

---

Published as a conference paper at ICLR 2026
The global loss function is defined as the cumulative mean squared error:
Protal
J
1o (T. Trarget dt
(E.2)
Here e00 is time indicator function. This function ensures that temperature errors computed only during the dish placement period.
E.1.2 GOVERNING EQUATIONS
The spatiotemporal evolution of the temperature field T(x, t) follows the two-dimensional inhomoge- neous heat conduction equation:
0T + XER?, 1
(E.3)
where a x 10 m the thermal diffusivity coefficient. The boundary condition adopts an infinite domain assumption: lim T(x, = To. The initial condition is set to the ambient temperature T(x, 0) To C The heat source term Q(x, () is constructed by superposing H, Gaussian- heating elements:
Q(x,)) t) = Σ
P(( j=1 PCp 2To²
H.
exp
x
12 20
(E.4)
Here a = 0.05 m is the effective radius of each heat source. The material properties are density p 7850 kg/m² and specific heat capacity Cp 460 J/(kg-K).
E.2 MULTI-RESOLUTION REINFORCEMENT LEARNING ENVIRONMENT
E.2.1
SPATIOTEMPORAL DISCRETIZATION SCHEME
We employ high-precision uniform Cartesian grid partition the physical domain Qt for bottom- layer physical evolution computation. The physical grid spacing is set to Arphy = 0.025 m. This corresponds Nphy x Nphhy 40 40 grid cells. The physical time step Atpky 1/450 < 0.0022 s. This satisfies the CFL stability condition for explicit finite difference schemes.
The coarsening degree of the RL observation space is controlled by parameter K 2 1. The RL grid size is N [Nphy/ K] with grid spacing Art = K Arphy- Spatial coarsening is implemented through regional averaging:
K-1K-1 ΣΣΣ p=0
T6:)
(Ks+p.K
(E.5)
K2
pky
9=0
where phy represents the temperature values on the physical grid. These values are obtained by evaluating the analytical solution at the corresponding spatial coordinates (i (iArphy p
The RL decision time step Al, is determined by parameter y = log (h,), where has Trotal/ Alphy is the temporal horizon at the physical layer and h, hy is the horizon at the RL layer. This design constructs a hicrarchical simulation structure. The underlying physical envi- ronment evolves H At Atp micro-steps for each decision executed by the agent. Different combinations of parameters K and y define 25 experimental configurations. The trade-off between sample efficiency and computational complexity across these configurations is the core focus of this study.
E.2.2 COTTE-ACTIO--RAEE
MECHANISM
The system state consists of the global temperature field T₁ E N.N.XNN, on the coarsened grid at the current time, along with the normalized remaining cooking time f, for each food item. Parameter K defines the coarse -rrsseiggainingg degree observations. Temperature readings are discretized at intervals of AT 2 x K. As K increases, the agent faces dual challenges of spatial information
43

---

Published as a conference paper at ICLR 2026
blurring (pixelation) and temperature quantization noise. These factors significantly intensify partial observability of the environment. The action space degrades correspondingly with the state space. Power adjustment of the two heat sources is modeled as discrete actions. The adjustment granularity directly coupled to through the minimum power increment AP = 50 W x K. This coupling causes the number of selectable power levels decrease significantly K increases. Coarse grids correspond to coarse control, which aligns with physical intuition.
To ensure consistency and physical meaning of the reward function across different time steps Arri we adopt continuous integration form based on physical micro-steps. Within the k-th RL decision step, the instantaneous reward rą is computed as the cumulative temperature error over all physical micro-steps in that period:
-
Numa Σ
(T.(kk,h) Tlarget To
Th Σ
10,,oook (bk,N)
Alphy
(E.6)
where H = Attt/Atphtt is the number of physical micro-steps within a single RL. step, and Tacale 50 the normalization constant. Through this approach, the cumulative reward R Σ remains an unbiased approximation of the original continuous control objective equation E.2 regardless of variations RL decision frequency. This design guarantees fairness in computational cost comparisons across different experimental groups.
E.3 TRAINING ALGORITHM AND CONVERGENCE CRITERIA
E.3.1 PPO ALGORITIIM IMPLEMENTATION
We adopt Proximal Policy Optimization (PPO) as the core training algorithm, considering the mixed characteristies of continuous state space and discrete action space in this problem. Compared to off-policy algorithms, PPO constrains policy update step sizes through clipped objective functions. This mechanism effectively prevents policy oscillations in complex non-stationary thermodynamic environments and ensures training robustness. Both the policy network (Actor) and value network (Critic) employ multi-layer perceptron (MLP) architectures with the global state vector as input. The hyperparameters are set as follows: learning rate 3 x 10 discount factor y = 0.99, and Generalized Advantage Estimation (GAE) parameter À = 0.95. Training proceeds synchronously across 8 parallel environments with GPU-accelerated gradient updates.
E.3.2 CONVERGENCE CRITERIA
To fairly compare sample efficiency across different resolution configurations (K, y). we adopt relative convergence criterion. Training considered converged and terminated when the agent achieves the preset performance threshold in three consecutive evaluations, where each evaluation based on the average over 10 complete episodes. The specific threshold is set theoretically derived optimal control baseline (see Section for the physical environment setup). This threshold setting ensures that computational cost comparisons across configurations are meaningful under comparable performance levels.
E.4 EXPERIMENTAL CONFIGURATION AND RESULTS
The experiment covers spatial resolution ratios K E {1.0, 2.0, 4.0, 6.0, 8.0) and temporal resolution parameter y log,.. (h,) E [1/8, 1/6, 1/4, 1/3, 1/2}. totaling 25 configuration points. Each con- figuration trained with five random seeds until reaching the convergence criterion. All experiments are executed on the same hardware platform (Intel Xeon Gold 6530, NVIDIA RTX 4090, 256 GB RAM).
E.4.1 S-N SCALING ANALYSIS
To validate the theoretical e-N scaling predictions, we analyze complete scaling curves for all 25 resolution configurations. Here, denotes the suboptimality gap defined e = R Reurents where R the optimal policy reward derived from the analytical solution. The variable represents training steps.
44

---

Published as a conference paper at ICLR 2026
Panel (c) of Figure shows E-N scaling curves for all 23 configurations that successfully reach the convergence threshold. Different colors represent temporal resolution y values, while color intensity indicates spatial resolution K (darker shades for smaller K, lighter for larger K). Shaded regions illustrate the performance range across different K values within each temporal resolution group. Vertical dashed lines mark the mean convergence step for each y, revealing clear separation between temporal resolution groups: finer temporal resolution (smaller y) requires fewer training steps to converge.
E.4.2 COMPUTATIONAL COST DECOMPOSITION
A key contribution of our experimental analysis is the rigorous decomposition of total computational cost into its constituent components. In AI4S reinforcement learning systems, computational re- sources consumed by two fundamentally different processes: (1) physical simulation on CPU, which involves solving or evaluating PDE-based environment dynamics, and (2) neural network operations GPU, which includes policy inference and gradient-based training updates.
This decomposition critical because CPU and GPU computational costs scale differently with resolution parameters. CPU costs are dominated by the complexity of the physical solver, which depends on grid resolution and time step constraints imposed by numerical stability. GPU costs are primarily determined by the number of RL training iterations required for convergence, which depends on the effective state-action space size and the difficulty of credit assignment under different temporal horizons.
Panels (d), (e), and (f) of Figure present heatmaps showing the CPU, GPU, and total computational costs respectively across all 25 configurations. Several important observations emerge from this decomposition:
CPU-GPU trade-off structure. The CPU and GPU cost landscapes exhibit distinct patterns. CPU costs tend to increase with finer temporal resolution (smaller loga... (h,)) )) due to the larger number of physical simulation steps required. GPU costs show more complex dependence on both spatial and temporal parameters, reflecting the interplay between state space complexity and horizon length in the RL training dynamics.
Non-additive composition. The total cost is not simply the sum of CPU and GPU costs in logarithmic scale. The relative contribution of each component varies significantly across the configuration space. At fine resolutions, CPU costs dominate due expensive physical simulations. At coarse resolutions, GPU costs become relatively more significant as the RI algorithm struggles with reduced observability and must compensate through additional training iterations.
Optimal configuration identification. The optimal configuration (K = 6.0. loga (h = 1/3)
achieves the minimum total cost by balancing these competing effects. This configuration provides sufficient resolution for the agent to learn effective control policies while avoiding the computational overhead of unnecessarily fine discretization.
E.4.3 COMPUTATIONAL COST SUMMARY
For each configuration, the logical interaction count Niotal recorded by environment destructors combined with the sampling multiplier Msampte (K, At) measured from post-training sampling experiments (see Section E.5). The equivalent total computational cost is reconstructed in units of
logio (a) Detailed numerical data with confidence intervals are presented in Table E.5
Table E.5: Equivalent Total Computational Cost (log to (FLOP Sweat )) with Confidence Intervals
logn.. (hr))
K =
K =
K = 4
K = 6
K 8
1/8
17.77 + 0.29
17.74 + 0.35
17.72 + 0.26
17.77 0.14
N/A N/A
1/6
17.59 + 0.42
0.29 17.60
17.65 + 0.21
17.65
1/4
17.59 + 0.26
17.56
17.50 + 0.17
17.54
17.81 0.27
1/3
17.64 +: 0.17
17.59 0.22
17.41 + 0.20
17.27+0.22
17.69 + 0.10
1/2
17.60 + 0.15
17.70
17.58 + 0.20
17.36 0.08
18.00 + 0.08
45

---

Published as a conference paper at ICLR 2026
Note: Values represent log1o ((oOP) with standard deviation computed from five random seeds. Both the logical interaction count Nostal during training and the statistical sampling multiplier Mample comprehensively considered. Configurations marked N/A fail to reach the performance threshold due to excessive discretization coarseness. Specifically, K' = 8.0 corresponds to only 5 x RL observation grids and action resolution of A P 400 W (with merely 3 selectable levels: W. 400 W, 800 W). Under this granularity, the agent cannot achieve fine temperature control. The optimal configuration (K 6.0, log (hr) 1/3) is highlighted in bold, corresponding
the lowest equivalent computational complexity of 10 0..2 FLOPs. This result confirms the theoretical prediction that both excessively coarse and excessively fine discretization lead to efficiency losses.
To ensure the reliability of our conclusions, we conducted extensive statistical analysis across multiple random seeds. Each configuration was trained with five independent random seeds, controlling for variations network initialization, environment stochasticity, and optimization trajectory. The confidence intervals reported in Table E.. and visualized in Figure 6 represent the standard deviation of the total computational cost across these runs.
E.5 COMPUTATIONAL COST RECONSTRUCTION METHODOLOGY
To ensure that experimental conclusions generalize to real A14S scenarios where numerical PDE solvers must be used, we adopt a computational cost reconstruction approach. The core challenge is that directly recording the CPU/GPU execution time for every state transition during training is complex and implementation--epennee computation (via destructor calls for object lifecycle management) may require vastly different convergence iterations depending on local condition numbers of the linear systems or nonlinear solver states at different grid cells. Recording these fine-grained execution times would require instrumenting every grid cell's 's solver, which is impractical for large-scale experiments.
Specifically, each invocation of the environment's physics
Instead, we adopt two-stage approach: first, we complete the full PPO training for all 25 resolution configurations and record only the logical interaction count Now via destructor counters in the environment core class. This count reflects the sample complexity of the RL algorithm itself, independent of the underlying physics solver implementation. Second, after training converges, we perform a post-hoc sampling study to estimate the computational multiplier that would be required in real scenarios using numerical solvers.
Specifically, we randomly select a representative subset of state-action pairs encountered during converged policy execution. For each selected state, we perform Monte Carlo sampling within the observation uncertainty range [x - Ay/2, x + Ay/2] to measure the frequency distribution of the next state falling into different grid cells. According to Theorem in the main text, the number of samples M Mamplee K2 A/) required to statistically distinguish the correct target grid with confidence
1 - 8 depends on the classification margin:
log(1/5) (8) ()lo((1/))
Msample ( K, At) =0
(E.7)
minn(Ap0)))
where Apli is the probability gap between the dominant grid cell and its closest competitor in dimension j.
We then estimate the per-sample computational cost (K))ver(K) for single PDE solver call at resolution K. For explicit finite difference methods with grid size Nphy/ K and CFL-constrained time step, this cost scales as Ssovver (K) ) ~ (Nphe//S) K)d per time step, where d is the spatial dimension. The final reconstructed total cost is:
FI.OPS K, At) = Niotal X Msampte K, At) x Cslvve(())
(E.8)
This reconstruction methodology ensures that reported computational costs honestly reflect what would be required in real A14S applications using numerical solvers, while still enabling efficient large-scale hyperparameter searches during the training phase. The sampling multiplier Msample accounts for the statistical overhead of state disambiguation under observation noise, and the solver computational burden. Similar reconstruction
cost Cselver accounts for the rewolution-deepeeeeeee can be applied to GPU costs by estimating the per-operation GPU FLOPs and memory bandwidth utilization, which typically differ from CPU costs by implementation speciffc constants but follow the same scaling laws.
46

---

Published as a conference paper at ICLR 2026
The physical environment observations in our experiments obtained from the analytical solution of the wwo-dimensional heat diffusion equation equation E.3 in free space. According to linear PDE theory, the temperature field evolution can be expressed through spatiotemporal convolution of the Green's function (heat kernel) with the heat source term:
T(x, t) = Tox + S/" S/g G(x,t; de dr
(E.9)
where the heat kernel C(x,, (; Ina(t-7] exp (- |x-E| unit point source in infinite space.
(-dx(1-7))
describes the diffusion pattern of
47


Reviews:
[{"review_id": "whF7eIIhBV", "reviewer": "Reviewer_QtUe", "summary": "The paper proposes a theoretical framework to determine the forward prediction error of RL systems control when trained with a neural surrogate. It also provides expressions for the optimal resolution ratio $K^*$ to minimize computational cost of different PDE systems. The authors verify their theoretical results by applying tabular RL to the Cart-Pole environment.", "questions": "- In the Cart-Pole experiment, do you train the RL with a neural surrogate model or just the physics simulator directly? Is this specified in the main text?\n- If one of the main problems with RL to control PDE based systems is the expensive cost of using simulators to train the RL, doesn't the same problem exist to train a surrogate neural emulator of the simulator? Are you assuming the amount of data required from a simulator for surrogate training less than that required for RL?\n- There are a few small misspellings in the paper which should be corrected:\n  - Line 357: \"To translate thees PAC-possible AI4S-RL analysis into practical system design guidance, this section establishes a computational resource allocation framework for AI4S-RL systems.\" Mispelling, should be \"To translate these\"\n  - Line 386, Theorem 3: \"yields optimal resolution raio between\" misspelled \"ratio\".", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The paper has some novelty when combining ideas from PAC learning, spectral analysis of PDEs, and coupling the sources of error from RL and AI4S systems.\n- The empirical results from the Cart-Pole experiment verify the results from their theoretical analysis.", "weaknesses": "- The paper only has empirical results for the Cart-Pole system, the theoretical results would be better justified if at least one other system was tested.\n- The paper assumes the only source of error from the neural surrogate is deterministic discritization error, but this is not always the case for real world systems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xDLtAwRIaY", "reviewer": "Reviewer_UFUv", "summary": "This paper addresses a fundamental and critical problem at the intersection of AI for Science (AI4S) and Reinforcement Learning (RL): how to manage the deterministic, resolution-dependent errors introduced when using AI-based surrogate models (like neural PDE solvers) as simulation environments for RL. Unlike stochastic noise, these discretization errors do not average out and can systematically bias the learned policy.\n\nThe authors' primary contribution is a novel theoretical framework, termed the \"$\\epsilon-N$ framework,\" designed to unify these two domains. This framework quantifies the minimal computational cost, $N^*(\\epsilon)$, required to train an RL agent to a value function accuracy of $\\epsilon$ with high probability.", "questions": "1. Could you please provide a more detailed justification for the $Cost(K)$ formulation in Theorem 3 and the \"computational balance condition\" $H_r^3 S_r A_r \\sim H_w S_w A_w$? An alternative might be $Cost_{Total} = N_{episodes}(\\rho) \\times H_r \\times Cost_{step}(K)$. How would the optimization of this (perhaps more direct) cost function compare to the results you presented?\n\n2. You empirically show that DQN benefits from a similar resolution balance. Theoretically, how do you hypothesize the $\\rho-K$ analysis would change when function approximation is introduced? Would the NN's approximation error add a new term, e.g., $\\rho = \\mathcal{O}(1/H_r + 1/K^d + \\epsilon_{NN})$? Or is the interaction more complex, perhaps with the NN's generalization ability reducing the required $S_r$ and thus changing the optimal $\\alpha$ and $\\beta$ in Theorem 3?\n\n3.  The analysis relies on a single $\\lambda_1$ for the system. In many of the systems studied (e.g., Tokamak, Airfoil), the instability $\\lambda_1$ is highly state-dependent (e.g., the plasma is more unstable in certain configurations). Does your framework assume a single \"worst-case\" $\\lambda_1$ for the entire state space? And if so, does this suggest that a truly optimal system would require adaptive resolutions ($K(s), H_r(s)$) that change based on the local, state-dependent $\\lambda_1(s)$?\n\n4. This framework's analysis is centered on discretization error ($\\rho$), which is assumed to vanish as resolution ($K, H_r$) increases. Nevertheless, the AI4S surrogates (such as neural operators), can also have systematic errors or model-form errors, if the surrogate is an imperfect representation of the true PDE over even infinite resolution (for example, due to a lack of training data or architectural limitations). How would this irreducible, non-vanishing error term modify the framework? Would it introduce a hard lower bound on the achievable accuracy $\\epsilon$, regardless of computational cost $N$?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The originality of this paper is outstanding. The problem itself is timely, but the authors' approach is, to my knowledge, new. Most work in this area either ignores surrogate error, treats it as generic bounded stochastic noise, or uses purely empirical tuning. This paper is the first I have seen to create a formal, unified theory that links PAC-RL sample complexity directly to the spectral and discretization properties of the underlying PDE. The $\\rho-K$ analysis and the resulting formula for $K^*$ in Theorem 3 are highly novel and insightful. The technical quality of the theoretical work is exceptionally high. The authors demonstrate a deep and fluent command of both reinforcement learning theory and advanced numerical analysis for PDEs. The error propagation analysis (e.g., in Section 4.2 and Appendix C.2.2) is incredibly detailed, correctly identifying non-obvious error sources like the sublinear $\\mathcal{O}(\\Delta a_w^{1/2})$ error from boundary control actuation and the reduced $H^{1/2}$ regularity at boundaries (via the trace theorem). Deriving concrete, closed-form scaling laws from this complex analysis is a major technical achievement.", "weaknesses": "* Focus on Tabular RL: The entire theoretical framework is built on tabular, finite-state-and-action PAC-MDP analysis ($S_r, A_r$). However, the motivating examples (Tokamak, Airfoil) are high-dimensional, continuous systems that would never be solvable with tabular methods; they inherently require deep function approximation (e.g., Deep Q-Learning). The paper provides an empirical result for DQN (Fig. 2) and notes it is \"consistent,\" but the theory does not extend to it. It is unclear how the Bellman error from function approximation would interact with the discretization errors ($\\rho$). This is a significant gap between the theory and its practical application.\n\n* Linear Dynamics Assumption ($\\lambda_1$): The framework's reliance on $\\lambda_1$, the dominant eigenvalue of the linearized operator, is elegant but is a simplification. Many complex systems (like turbulence in the airfoil example) are dominated by nonlinear error growth, not just linear modal instability.\n\n* Vagueness in Cost Function Definition: The definition of the final $Cost(K)$ in Theorem 3 is slightly confusing. It appears to be $Cost_{RL} \\times (\\text{Penalty for } \\rho)^2$, where $Cost_{RL}$ is the tabular RL sample complexity $\\mathcal{O}(H_r^3 S_r A_r)$. This $Cost_{RL}$ term already scales with $K$ (via $S_r \\sim K^\\alpha, A_r \\sim K^\\beta$). However, the cost of a single surrogate call (which should be the primary driver of $N$ and scale with the AI4S grid, i.e., $\\mathcal{O}(K^d)$) seems absent. The paper instead posits a \"computational balance condition\" $H_r^3 S_r A_r \\sim H_w S_w A_w$ to relate the RL cost to the (absent) simulator cost. This feels like an assumption to make the optimization tractable, rather than a derivation of total computational cost.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1aNwUqCbKk", "reviewer": "Reviewer_zDg6", "summary": "This work is concerned with analyzing the error introduced by using approximate simulators being used in solving complex physical control problems, governed by PDEs. The authors propose a principled framework, namely $\\epsilon-N$ framework to quantify the minimum computational cost required for such approximate simulators to ensure that tabular RL estimates the value function without any bias, with a large probability. Through their theoretical framework, the authors provide optimal resolution trade-offs for some of the AI4S-RL systems such as Tokamak, and Cart-Pole.", "questions": "1. I understand that the discretization scheme we use will impact the solution we obtain. But, if we use continuous state and time representations, and assuming $\\Delta t \\approx 0$ (for state evolution),  do we face the similar issues? At least in the case of Cart-Pole, $\\Delta t \\approx 0$ doesn’t pose a huge computational burden. Is that correct?\n\n2. Perhaps just use $w$ to denote the AI4S domain? $x_p$ just appears for one case; might be easier to read.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. It is important to understand how accurate the RL environment should be to compute optimal policy for the actual physical system, and this work tries to answer that. \n\n2. The paper is relatively easy to follow despite being heavy on theory.", "weaknesses": "1. I would have liked to see the empirical result in the main paper instead of it being in the appendix. \n\n2. I think a careful rearrangement of symbols is necessary as some of the symbols just appear without any definition. For example, $\\lambda_1$ in line 374 is mentioned without context, although it is later mentioned in the appendix. $\\lambda$, being a crucial parameter, might benefit from a short introduction.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Vqb0LlWee5", "reviewer": "Reviewer_ktCj", "summary": "The paper proposes a theoretical framework connecting numerical resolution (in AI4Science simulators) and decision resolution (in RL agents). It claims an analytical form for the optimal resolution ratio and derives scaling laws linking discretization errors to computational cost. Experiments on a learned Cart-Pole dynamics model show non-monotonic performance trends with respect to resolution, suggesting the existence of an optimal middle point.", "questions": "How was the spectral constant estimated in practice? \n\nCan you provide at least one full experimental validation on a real PDE-based simulator (e.g., tokamak or CFD)?\n\nHow sensitive are the results to random seeds and network initialization?\n\nDoes the proposed cost function (Eq. 12) correlate with real compute cost (GPU-hours, memory, etc.)?\n\nCould the N^{1.6} exponent be confirmed by sampling multiple off-optimal configurations instead of comparing two extremes?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper bridges numerical simulation fidelity and RL resolution, a connection that is conceptually valuable for AI4Science.\n\nThe Cart-Pole experiments show that increasing resolution indefinitely does not always improve performance; there exists an intermediate optimum.", "weaknesses": "The central results (universal scaling, closed-form K* and N^1.6 deviation cost) are not empirically verified beyond the Cart-Pole example.\n\nThe paper lists tokamak, airfoil, and teppanyaki as examples but provides only theoretical discussion, no numerical results.\n\nThe Cart-Pole environment itself is a learned world model, so conclusions about “simulator–RL balance” may not directly generalize to physical simulators.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)