Distance-Residual Physics-Informed Neural Networks: A Deep Learning Framework for Differential and Partial Differential Inclusions
Authors: Maria Filipkovska, Juan José Marín, Isil Oner, Francisco Periago, Krzysztof Rykaczewski
Organizations: B. Verkin Institute for Low Temperature Physics and Engineering of the National Academy of Sciences of Ukraine, Nauky Ave. 47, Kharkiv 61103, Ukraine. · Department of Applied Mathematics and Statistics, Universidad Politécnica de Cartagena, member of the European University of Technology EUT+, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain. · Department of Mathematics, Faculty of Sciences, Gebze Technical University, Kocaeli 41400, Turkey. · Faculty of Mathematics and Computer Science, Nicolaus Copernicus University in Toruń, Chopina 12/18, 87-100 Toruń, Poland.
We introduce Distance-Residual Physics-Informed Neural Networks (DR-PINNs), a physics-informed learning framework for approximating solutions of ordinary and partial differential inclusions (DIs), governing laws in which a differential operator is constrained to lie in a set-valued map rather than equaling a prescribed function. The method replaces the classical pointwise PDE/ODE residual by the squared distance from the differential operator to the admissible set. This distance vanishes exactly when the inclusion is satisfied and measures the infimal correction needed for the operator to enter the admissible set. For a fixed closed convex admissible set, the squared distance is differentiable with respect to the operator value, with gradient given by the metric projection. When the admissible set also depends on the network state, that dependence is included through the chain rule. The framework encompasses ordinary and partial DIs with set-valued reaction terms. For both settings we prove consistency: under the stated closedness, measurability, convexity, and growth assumptions, any sequence of candidates satisfying the initial (and, in the parabolic case, boundary) conditions whose continuous distance-residual functional tends to zero admits a subsequence converging to an exact solution of the target inclusion. These are conditional statements for the continuous distance-residual functional; they do not cover the finite-collocation training loss, the behavior of the optimizer, or convergence rates. For admissible sets given as convex hulls of finitely many vertices, projection onto the set reduces to a small convex quadratic program, making the loss efficiently computable inside the training loop. Numerical experiments demonstrate high accuracy on the differential-inclusion benchmarks.
Figures & tables
Figure 1: Raw vertex set Uraw ( Kraw=12 , blue circles: the four true extreme vertices of the ℓ1 diamond plus eight redundant interior points) and the four extreme vertices K=4 recovered by Quickhull (red stars), together with the convex hull of Uraw .
Figure 2: Training history of the inclusion loss L(θ) for the single DR-PINN of this section (semi-logarithmic scale).
Figure 3: Phase portrait (left) and state components versus time (right) for the trained trajectory xθ(t) , compared with the sampled reference cloud Rref(t) (304 trajectories, light blue) and one representative bang-bang extremal (dark blue).
Kraw
8
16
32
64
128
K
6
6
6
6
6
Table 1: Extreme vertices K recovered by Quickhull as Kraw grows, at fixed d=2 .
d
2
4
6
8
K
12
23
32
32
Table 2: Extreme vertices K recovered by Quickhull as the state dimension d grows, at fixed Kraw=32 .
Figure 4: Per-epoch training time with (blue) and without (red) the offline Quickhull reduction; markers are means over 5 repetitions after a discarded warm-up run, error bars show one standard deviation. Left: as a function of Kraw at fixed d=2 (Table 1 ). Right: as a function of the state dimension d at fixed Kraw=32 (Table 2 ).
λξ
∣xz(T)−xf∣
mean ℓ
max ℓ
10−3
2.0×10−5
0.972
0.982
10−2
4.8×10−6
0.971
0.982
10−1
9.7×10−6
0.968
0.982
Table 3: Continuation over the selector-regularization weight λξ (warm-started L-BFGS-B): terminal mismatch and level statistics of the optimized selector. Admissibility ℓ(ti,ξz(ti))<1 holds at every grid point and every stage, by construction.
Figure 5: Admissible velocity tube of the inclusion ( 6.9 ) in (v1,v2,t) space. The surface is the boundary of the translated rotating ellipses g(t,xz(t))+E(t) ; the blue curve is the centerline g(t,xz(t)) , the orange curve the selected velocity x˙z(t) , and the grey segments the selector displacements ξz(ti) .
Figure 6: Optimized selector ξz(t) inside the rotating control ellipse E(t) (thin blue rings, drawn every sixth grid point). The color of the points encodes the level ℓ(ti,ξz(ti)) , which remains below 1 at the plotted time-grid nodes by construction.
Figure 7: Left: optimized state trajectory of ( 6.11 ) from x0=(0,0)⊤ to the target xf=(1.8,0.1)⊤ (final continuation stage λξ=10−1 , terminal mismatch 9.7×10−6 ). Right: admissibility level ℓ(t,ξz(t)) along the horizon; the dashed line marks the boundary value ℓ=1 of ∂E(t) .
λ
grid
Δt
t∗
Δt∗
max∣Δ∥u∥∣
max∣Δupre∣
0.4
49×49
0.0005
0.1845
+0.0000
0
0
0.4
49×49
0.00025
0.1840
-0.0005
4.5×10−4
1.7×10−4
0.4
49×49
0.000125
0.1839
-0.0006
6.7×10−4
1.1×10−4
0.4
97×97
0.0005
0.1845
+0.0000
5.1×10−5
4.1×10−5
0.4
97×97
0.00025
0.1840
-0.0005
5.0×10−4
1.3×10−4
0.4
193×193
0.00025
0.1840
-0.0005
5.1×10−4
1.2×10−4
Table 4: Refinement study for the reference solver: detected extinction time t∗ (first grid time with an exactly zero discrete field; persistent up to T in all runs), its change Δt∗ with respect to the baseline configuration ( 49×49 , Δt=5×10−4 ), the maximum change of the L2(Ω) -norm decay curve, and the maximum change of the solution profile shortly before extinction (at 0.9t∗ , interpolated to the baseline grid). At fixed Δt the detected t∗ is independent of the mesh; the variation under Δt -refinement is of the order of one baseline time step and is dominated by the temporal quantization of the detection.
Figure 8: Reference solution: ∥u(t)∥L2(Ω) on a log scale for λ∈{0,0.4,0.6,0.8} . The vertical dotted lines mark the finite-time extinction instants t∗ of the three relay cases; the λ=0 heat-equation curve does not extinguish within [0,T] . Curves are clipped below 10−6 for readability; after extinction the reference norm is exactly zero (drop below the axis).
Figure 9: Training history of the inclusion loss for all 45 runs of the campaign, colored by ε (validation on a fixed batch of 4000 points). All banded runs ( ε>0 ) escape the initial plateau and continue to decrease; the five exact-relay runs ( ε=0 , top curves) stagnate an order of magnitude higher.
λ
tε,ref
tε (DR-PINN) median [range]
persistent crossings
violations (of 4000 )
tail dist. to 0
witness admissible
0.4
0.1350
0.1350 [0.1350, 0.1375]
5/5
920
1.8×10−3
yes
0.6
0.1275
0.1275 [0.1250, 0.1275]
5/5
896
1.8×10−3
yes
0.8
0.1200
0.1200 [0.1200, 0.1225]
5/5
867
4.0×10−3
no
Table 5: λ -sweep at ε=0.05 (front sampler, 5 seeds per row; violations and tail distances are seed medians). The DR-PINN band-entry time ( 6.28 ) agrees with the reference band-entry time to within one evaluation-grid step for every seed, every first crossing is persistent, and the monotone ordering in λ is preserved by every individual run — in contrast to the raw 10−3 threshold crossings of the exact relay reported as the ε=0 control in Table 6 . “Witness admissible” records whether λ∥φ∥L∞≤ε (Remark 6.3 ); in all cases the trained tail is an order of magnitude closer to the zero state than to the extremal witness (cf. Table 6 ), so the networks select the extinct solution even where the ε -banded funnel is nontrivial.
ε
tε,ref
tε (DR-PINN) median [range]
persist. cross.
final val. mean
violations (of 4000 )
tail dist. to 0
tail dist. to wλ
0
—
— ( 2/5 no crossing)
0/5
0.437
3 911
1.2×10−3
2.6×10−2
0.01
0.1550
0.1550 [0.1475, 0.1575]
4/5†
0.081
1 301
2.4×10−3
2.4×10−2
0.02
0.1475
0.1475 [0.1450, 0.1475]
5/5
0.047
1 152
1.8×10−3
2.3×10−2
0.05
0.1275
0.1275 [0.1250, 0.1275]
5/5
0.027
896
1.8×10−3
2.5×10−2
0.1
0.1050
0.1050 [0.1050, 0.1050]
5/5
0.017
652
2.9×10−3
2.5×10−2
Table 6: ε -sweep at λ=0.6 (front sampler, 5 seeds per row; final validation mean, violations, and tail distances are seed medians). The exact relay ε=0 is the negative control and reproduces the failure mode of the discontinuous target: no persistent crossings, an order-of-magnitude larger residual, and near-total violation of the pointwise inclusion on the evaluation cloud. From ε=0.02 on, band-entry times match the reference to within one grid step with persistent crossings for all seeds. † At ε=0.01 one seed’s first crossing is transient; its persistent crossing occurs only at t=0.25 . The tail is consistently an order of magnitude closer to u≡0 than to the extremal witness wλ of Remark 6.3 .
Figure 10: Band-entry times across λ∈{0.4,0.6,0.8} at ε=0.05 : DR-PINN values (five seeds per λ , dots) against the reference band-entry times and the reference extinction times t∗ (stars). The DR-PINN values sit on the reference band-entry values to within the evaluation-grid resolution and preserve the monotone ordering in λ run by run.
Figure 11: Dense-grid distance-residual statistics along the ε -sweep at λ=0.6 ( ε=0 plotted at 2×10−3 on the logarithmic axis). The mean residual decreases monotonically with ε and drops below the λ2 plateau of the smooth-ansatz bound as soon as the band opens, while the maximum residual — attained ahead of the moving interface — remains of order λ2 times the squared operator overshoot for every ε .
Figure 12: L2(Ω) -norm decay for λ=0.6 , ε=0.05 , on a log scale: reference solution (solid) and the five DR-PINN seeds (dashed). The horizontal dash-dotted line marks the band level ε , the vertical dotted line the reference extinction time tref∗=0.1645 . The curves track the reference to within ≈2% before band entry; after band entry the smooth surrogates retain a small nonzero tail, which the band-entry observable is designed to ignore.
Figure 13: Branch-selection diagnostic ( λ=0.6 , ε=0.05 ): empirical cloud of (uθ,∂tuθ−Δuθ) at 4000 post-training collocation points, overlaid on the graph of Φε (red; the shaded step of width 2ε is the band). Points with s>ε cluster tightly on the branch z=−0.6 ; inside the band the residual vanishes; the vertical scatter concentrates in a thin layer ahead of the moving interface, where the discontinuity of the relay limits pointwise accuracy. The inset magnifies the band region.
sampler
tε median [range]
persistent crossings
final val. mean median [range]
violations median
front (adaptive)
0.1275 [0.1250, 0.1275]
5/5
0.027 [0.020, 0.041]
896
RAR
0.1275 [0.1275, 0.1275]
5/5
0.023 [0.017, 0.026]
899
uniform
0.1275 [0.1250, 0.1300]
5/5
0.026 [0.016, 0.027]
895
Table 7: Sampler ablation at (λ,ε)=(0.6,0.05) , 5 seeds per sampler. The three schemes are statistically indistinguishable on every reported metric; the band-entry observable is insensitive to the collocation strategy at these settings.
Physics-Informed Neural Networks (PINNs) recast PDE solving as an optimisation problem in function space by minimising a residual-based objective, yet many applications require additional derivative-based relations that are just as fundamental as the governing equations. In this paper, we present Derivative-Constrained PINNs (DC-PINNs), a general framework that treats constrained PDE solving as an optimisation guided by a minimum objective function criterion where the physics resides in the minimum principle. DC-PINNs embed general nonlinear constraints on states and derivatives, e.g., bounds, monotonicity, convexity, incompressibility, computed efficiently via automatic differentiation, and they employ self-adaptive loss balancing to tune the influence of each objective, reducing reliance on manual hyperparameters and problem-specific architectures. DC-PINNs consistently reduce constraint violations and improve physical fidelity versus baseline PINN variants, representative hard-constraint formulations on benchmarks, including heat diffusion with bounds, financial volatilities with arbitrage-free, and fluid flow with vortices shed. Explicitly encoding derivative constraints stabilises training and steers optimisation toward physically admissible minima even when the PDE residual alone is small, providing reliable solutions of constrained PDEs grounded in energy minimum principles.
Kentaro Hoshisashi, Carolyn E Phelan, Paolo Barucca
Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK
Physics-Informed Neural Networks (PINNs) are a machine learning method for solving forward and inverse Partial Differential Equations (PDEs). When applied to PDEs with Dirac delta functions in the forcing terms, boundary conditions, or initial conditions, PINNs require approximating them with smooth surrogate functions, a practice that can introduce significant modeling errors. In this work, we exploit the interpretation of PINNs as Residual Least Squares (RLS) methods and show that this perspective enables direct treatment of Dirac delta terms by integrating the weak-form equation. Among RLS formulations other than PINN, we focus on the Radial Basis Function (RBF) expansion (also known as a single-layer RBF Network). We show that while integrating out the Dirac delta in PINNs causes residuals to fail to converge to zero, RBF-RLS consistently provides good forward and inverse solutions to transport problems. We explain this finding using the Neural Tangent Kernel (NTK) theory. We test both approaches on linear PDEs that represent groundwater flow and transport in porous media and rivers. We solve inverse problems to fit synthetic data, noisy synthetic data, and real-world measurements.
Manuel Reyna, Alexandre Tartakovsky
Department of Civil and Environmental Engineering, University of Illinois Urbana-Champaign, 205 North Mathews Avenue, Urbana, 61801, IL, USA
Neural operators evaluate parametric partial differential equations cheaply but degrade sharply outside their training distribution. Physics-informed neural networks avoid dependence on labeled data, yet their optimization can be basin-fragile: when the governing residual admits multiple solutions, a PINN trained from scratch may converge to a physically incorrect state despite achieving a small residual. We show that these failure modes can be addressed jointly: an imperfect NO provides the structural prior needed to place a PINN in the correct solution basin, while the PDE residual refines the solution beyond the operator's accuracy. We introduce a three-stage framework that freezes the spatial basis of a physics-informed NO, extrapolates its solution branch to an out-of-distribution parameter using a polynomial continuation prior, and distills the resulting field into a fresh PINN. The NO need not be accurate at the target; it transfers solution-branch information, while PDE residual minimization in the PINN governs convergence. We evaluate the framework on three nonlinear PDEs: 1D viscous Burgers, 2D steady Allen-Cahn near a pitchfork bifurcation, and 2D steady lid-driven cavity flow. For Allen-Cahn, where the trivial solution satisfies the PDE residual exactly, a standard PINN collapses to the trivial zero branch, whereas distillation from the crude extrapolated operator recovers the non-trivial branch that matches the finite-difference reference. For the lid-driven cavity, extrapolating to a Reynolds number of Re = 3200 accelerates convergence to the correct physical state, achieving competitive accuracy using fewer parameters and optimization steps than recent literature baselines. These results establish a simple principle: an NO need not accurately predict the solution to be useful; it only needs to identify the correct basin from which PINN optimization can recover it.
S. Mohammad Mousavi, Teeratorn Kadeethum, Nikolaos Bouklas +1
Sibley School of Mechanical and Aerospace Engineering, Cornell University · AI Lab, Siemens Energy · Department of Civil and Systems Engineering, Johns Hopkins University