Full text
SCAN 2025 20th International Symposium on Scientific Computing, Computer Arithmetic, and Verified Numerical Computations Book of Abstracts Editors: Ekaterina Auer1 Marit Lahme2 Andreas Rauh2 1University of Applied Sciences Wismar, Department of Electrical Engineering and Computer Science, Wismar, Germany 2Carl von Ossietzky Universität Oldenburg Department of Computer Science Distributed Control in Interconnected Systems, Oldenburg, Germany Oldenburg, Germany, September 22-26, 2025 Universität Oldenburg / Matthias Hornung
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany I Aim of SCAN 2025 The series of International Symposia on Scientific Computing, Computer Arithmetic, and Verified Numerical Computations (SCAN) has been continued with the 20th edition from the 22nd to 26th of September 2025 in Germany. We have been pleased to invite around 50 international participants from the Czech Republic, France, Germany, Hungary, Japan, Poland, Switzerland, and the United States to the Carl von Ossietzky Universit¨at in the city of Oldenburg, a distinguished location in the German state of Lower Saxony. SCAN 2025 has been an excellent meeting place for researchers from the fields of reliable computing, software engineering, and uncertainty quantification and those from such wide and varied areas as robotics, control, structural and civil engineering, and signal processing. Here, long-standing colleagues have been able to meet again in a relaxed though productive setting after a long absence of an in-person SCAN event. The new edition of the conference has continued to further strengthen the exchange of novel scientific ideas and contained not only classical presentations in a lecture format but also discussion sessions focused on young researchers. In those, PhD students have has the opportunity to present their research activities in more detail and exchanged views about them with a broad audience of more experienced scientists or other PhD students. The scientific committee of SCAN has decided in its meeting that the next conference (SCAN 2027) will be held in Kitakyushu, Japan; SCAN 2029 is planned to be held in Krakow, Poland.
Contents Monday, September 22, 2025 – Library Auditorium 1 Constructing the Bessel function rigorously via the power series arithmetic H. Miyauchi, T. Asai, M. Kashiwagi, and A. Takayasu 2 Interval Uniform, Non-Uniform, Rational, Non-Rational B-spline Curves L. Si Larbi, E. Lucet, and J. Alexandre dit Sandretto 5 Tuesday, September 23, 2025 – Library Auditorium 8 Verified Error Bounds for Sparse Systems S.M. Rump 9 Formal Verification of State and Temporal Properties of Neural NetworkControlled Systems A. Besset, J. Tillet, and J. Alexandre dit Sandretto 10 Exploiting the Impossible: Towards Resilience of Decision Making Against Misperceptions M. Fr¨anzle, P. Kr¨oger, and A. Nienaber 13 Computer-assisted proof of the simplicity of the second Dirichlet eigenvalue for non-equilateral triangles R. Endo and X. Liu 15 Tuesday, September 23, 2025 – A01-0-006 17 Verry: an open-source package for verified computation written in Python 3 R. Iwanami 18 Hardware accelerated interval arithmetic for mobile robotics using RISC-V ISA extension P. Filiol, T. Bollengier, L. Jaulin, and J.-C. Le Lann 20 Towards Interval Arithmetic in TensorFlow: A Comparison of Approaches J. Khun and J. Schmidt 22 Wednesday, September 24, 2025 – Library Auditorium 25 Interval and Set-Based Approaches for Control and State Estimation: Their Use for Offline and Online Purposes A. Rauh 26 II
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany III Green-Representable Solutions: Reformulating Suband Super-solution Theory for Poisson’s Equation K. Tanaka, R. Iwanami, K. Matsue, and H. Ochiai 29 Verified Computation of All Positive Solutions to a H´enon-Type Boundary Value Problem T. Asai, K. Tanaka, S. Tanaka, and S. Oishi 31 Semigroup approach for validating solutions to semilinear parabolic PDEs A. Takayasu and J.-P. Lessard 33 Wednesday, September 24, 2025 – A01-0-006 35 For statistical analysis of big data, interval uncertainty is needed O. Kosheleva and V. Kreinovich 36 How to compare situations in which we measure different quantities with different uncertainty J. Alam, I. Hossain, T. Hossain, M.N. Sojib, O. Kosheleva, and V. Kreinovich 38 Towards Fair and Explainable Medical Risk Prediction Software via Dempster-Shafer Theory E. Auer and W. Luther 40 Thursday, September 25, 2025 – Library Auditorium 42 Bringing Formal Methods from Academia to Real-World Applications in Industry – My Personal Two-Decades-Journey – T. Teige 43 Separator for the remoteness constraint Q. Brateau, F. Le Bars, and L. Jaulin 44 Computing Interval Detection-Probability Grids via Monte-Carlo method for Underwater Robotics D. Esnault, S. Rohou, F. Le Bars, and L. Jaulin 47 Set-Based Identification of Characteristic Curves and Its Challenges for Real-Life Applications M. Lahme and A. Rauh 49 Fast Implementation of Interval Matrix Multiplication Using InfimumSupremum Representation With SIMD Operations H. Kijima and T. Ogita 51 Verification of Singular Values under Oblique Inner Product Space for Matrices T. Terao, Y. Watanabe, and K. Ozaki 53
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany IV Tight Enclosure of Matrix Multiplication using Fused Multiply-Add K. Ozaki and T. Koizumi 55 Parameter Robustness of Neural Networks A. Sz´asz and, B. B´anhelyi 57 Interval Based Verification of Adversarial Example Free Zones for Neural Networks T. Csendes 59 Thursday, September 25, 2025 – A01-0-006 61 Inconsistencies in Fuzzy Estimations: Kaucher Arithmetic Naturally Appears O. Kosheleva and V. Kreinovich 62 Shapley Value Under Interval Uncertainty Revisited: Why Seemingly Natural Axiomatic Approach Is Not Fully Adequate M. Svitek, O. Kosheleva, and V. Kreinovich 64 A new wrapper for a reliable resolution of underdetermined nonlinear equations L. Jaulin 66 Linear Programming Problems with Absolute Values and Interval Uncertainty M. Hlad´ık 68 Basis stability in interval quadratic programs C. Koteck´y and M. Hlad´ık 70 B&P algorithms for continuous constraint problems: a survey of branching strategies C. Jermann, N. Revol, and C. Solnon 72 Estimation of the Domain of Attraction for Nonlinear Systems using the Bihari Inequality R. Dilji, B. Tibken, R. Dehnert, Y. Fan, R. Deisling, and L. Ackerschott 74 Observer-Based Approaches for a Verified Simulation and Pseudo State Estimation of Fractional Dynamic Systems A. Rauh and M. Lahme 77 Friday, September 26, 2025 – Library Auditorium 79 Automated Verification of Discrete Probabilistic Programs C. Matheja 80
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany V Investigating chaos in Delay Differential Equations with rigorous numerical methods A. Gierzkiewicz, J. Kural, and R. Szczelina 82 Computer assisted proof of existence of periodic solutions to ENSO delay differential equation model J. Kural, A. Gierzkiewicz, and R. Szczelina 84 Ultra-wideband Based Smart Wheelchair Static Pose Estimation using Interval Analysis T. Le Terrier, M. Babel, and V. Drevelle 86 Set-Based Contracts for Systematic Controller Tuning in Interconnected Dynamic Systems A. Rauh and F. Bruns 89 Oscillating orbits in the Sitnikov model: equal masses case M.J. Capi´nski, A. Gierzkiewicz, and P. Mart´ın 91 Interval Particle Filter for LiDAR-Based Object Tracking M. Fnadi and R. Lherbier 93 Friday, September 26, 2025 – A01-0-006 95 Improvements of the Geometrical Test in Interval Branch and Bound methods M. Gencsi, B. G.-T´oth 96 Exhaustive Interval-based 2-D Shape Registration Under Similarity Transformation V. Radwan, S. Rohou, and G. Trombettoni 98 Adaptative parallelepipedic approximation of the image of a set by a nonlinear function M. Godard, L. Jaulin, and D. Mass´e 100 High-Performance Emulation of Matrix Multiplication using INT8 Matrix Engines and its Error Analysis Y. Uchino, K. Ozaki, and T. Imamura 103 Efficient Acceleration Strategies for Interval Branch-and-Bound Type Methods L. Gillner and E. Auer 105 GPU-Accelerated Algorithmic Differentiation For Reliable Computing: Comparing Different Architectures D. Romano, E. Auer, F. Gregoretti, and L. Gillner 108
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 1 Monday, September 22, 2025 Library Auditorium 8:30–12:20 CoProD: 17th International Workshop on Constraint Programming and Decision Making 12:30–14:00 Lunch Break – Food truck 14:00–14:30 Opening 14:30–16:00 Moore Prize Lecture Tristan Buckmaster, Gonzalo Cao-Labora, Javier Gomez-Serrano: Smooth imploding solutions for 3D compressible fluids 16:00–16:30 Coffee Break 16:30–17:00 Regular Session: Special Functions 16:30–17:00 Hiroaki Miyauchi, Taisei Asai, Masahide Kashiwagi: and Akitoshi Takayasu: Constructing the Bessel function rigorously via the power series arithmetic 17:00–17:30 Lucas Si Larbi, Eric Lucet and Julien Alexandre Dit Sandretto: Interval Uniform, Non-Uniform, Rational, Non-Rational Bspline Curves 17:30–18:30 PhD Poster Session All posters presented in this PhD session are included subsequently together with the associated abstract.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 2 Constructing the Bessel function rigorously via the power series arithmetic Hiroaki Miyauchi1, Taisei Asai2, Masahide Kashiwagi2, Akitoshi Takayasu3 1Graduate School of Science and Technology, University of Tsukuba [email protected] 2Faculty of Science and Engineering, Waseda University [email protected] [email protected] 3Institute of Systems and Information Engineering, University of Tsukuba [email protected] Keywords: Power series arithmetic, Bessel function, Interval arithmetic Introduction This study introduces a verified numerics framework for integrals involving the Bessel function that arise in the analysis of nonlinear elliptic boundary value problems. Existing quadrature schemes rarely control rounding and truncation errors rigorously, and the accuracy of their results cannot be guaranteed. In particular, to the best of our knowledge no existing approach encloses in interval form integrals of the type Z1 0 p+1 Y i=1 Jni(νni,mir)rdr, ni= 0,1,2, . . . , mi= 1,2,3,..., where Jn(x) is the n-th order Bessel function of the first kind and νn,m denotes the m-th positive root of Jn(x). This integral corresponds to the Galerkin projection in semilinear elliptic equations (1) shown at the end of this abstract. There exist methods for rigorously computing the Bessel function, for example, within the Arb library [1]. This is a C library for rigorous real and complex arithmetic with arbitrary precision based on ball arithmetic that implements algorithms for computing Jn(x) and Yn(x) (the n-th order Bessel function of the second kind) in the real and complex domains. However, Arb lacks built-in tools for verifying quadratures of the Bessel function. This makes it difficult to account for truncation errors in such quadratures directly. With the background mentioned above, we have developed a method that performs verified numerics for integrals involving the Bessel function, enclosing the integral values directly and rigorously in interval form. In this talk, providing the Bessel function using the power series arithmetic in the kv library [2], a collection of C++ libraries for verified numerical computations, we can obtain rigorous enclosures of the Bessel function via interval arithmetic. We apply the Schl¨omilch’s diffrential recurrence formula [3] to construct the power series expansion of Jn(x). All coefficients of the power series up to degree nare rigorously included in the interval coefficients. The remainder term is also rigorously included
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 3 in the coefficient of the highest degree. This yields a verified power series representation of Jn(x). Additionally, the verified power series representation provides interval inclusion of the values of integrals such as Z1 0 J0(ν0,1r)rdr, Z1 0 J0(ν0,1r)3rdr. Selected results in practice Furthermore, using the power series expression, the verified quadrature of the Bessel function is rigorously included as follows: Z1 0 J0(ν0,1r)rdr ∈h0.21587740350984231 3808i Z1 0 J0(ν0,1r)3rdr ∈h0.097461301068597087 4366i Z1 0 J1(ν0,1r)J1(ν0,2r)rdr ∈h−2.35055030994e−15,2.39776731456e−15i Z1 0 J0(ν0,1r)J1(ν1,1r)J2(ν2,1r)rdr ∈h0.036053739292458323 183674i Z1 0 J8(ν8,1r)J9(ν9,1r)J10(ν10,1r)rdr ∈h0.0046401195373457745 52999569i, where the interval [0,1] is partitioned into sixteen sub-intervals and the Bessel function is expanded by the power series up to the degree 20. For applications of the provided method, we aim to consider verified numerics for solutions to semilinear elliptic equations on the unit disk. −∆u=f(u) in Ω u= 0 on ∂Ω,(1) where Ω = {(r, θ)|0≤r≤1,0≤θ≤2π} ∈ R2,f(u) is a polynomial of order p. Since the Bessel function is the eigenfunction of the Laplacian on the unit disk, we can expect a highly accurate approximation of the boundary value problems. To this end, the methods for integrals introduced in this abstract are essential. This prospective application to the problem in (1) serves as the primary motivation for this study. References [1] Johansson, F: Arb kibrary (2013). https://arblib.org/index.html [2] Kashiwagi, M.: kv library (2016). http://verifiedby.me/kv/index-e.html [3] Watson, G.N: A Treatise on The Theory of Bessel Functions, Cambridge (1966).
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 10 Formal Verification of State and Temporal Properties of Neural Network-Controlled Systems Antoine Besset1, Joris Tillet1and Julien Alexandre dit Sandretto1 1ENSTA Paris, Institut Polytechnique U2IS Palaiseau, France {antoine.besset,joris.tillet,alexandre}@ensta.fr Keywords: Signal Temporal Logic, Verification of Neural Network, Interval Methods, Cyber-Physical Systems Ensuring the safety of Neural Network Controlled Systems (NNCS) remains a major challenge due to the opaque nature of neural networks, especially when temporal properties are involved. This paper presents an interval analysis-based framework for verifying both state and temporal properties of NNCS using Signal Temporal Logic (STL) specifications [?,?,?,?]. We introduce an STL monitoring algorithm based on interval analysis, featuring adaptive time sampling and formal guarantees of satisfaction over continuous domains. The STL formalism allows a rich temporal property specification while our approach is broadly applicable to neural networks when activation functions can be expressed as Ordinary or Differential Algebraic Equations (ODEs/DAEs). Reachability analysis, following the differential approach of [?], employs an ODE solver with affine arithmetic to ensure tight enclosures and dependency tracking. We demonstrate the effectiveness of the method on two case studies, involving both NNCS and systems with temporal constraints. To express temporal properties, we adopt a temporal logic formalism known as Signal Temporal Logic (STL) [?,?]. It has been applied in the domains of robotics and control. STL formulas allow the expression of various temporal properties using explicit time bounds, combining logical connectives with bounded Until temporal operators (U[a,b]) [?]. The syntax of STL is defined recursively as follows: ϕ:= µ|T| ¬ϕ|ϕ1∨ϕ2|ϕ1U[a,b]ϕ2.(1) We extend the verification of predicate (µ) with an inclusion predicate (Xµ) to verify properties on reachable tubes y(t, [y0]) ⊆([˜y], t),∀t∈[t0, T],[?,?]. A reachable set at tis [˜y](t). ([˜y], t)⊨µi:= 1,if [˜y](t)⊂ Xµ, 0,if [˜y](t)∩Xµ=∅, [0,1],otherwise. (2) To evaluate satisfaction, we extend the STL syntax with the Boolean Interval Arithmetic [?,?], enabling sound reasoning under uncertainty. To conduct reachability analysis of an NNCS, one effective approach is to exploit the differential properties of its activation functions [?]. For instance, in the case of the sigmoid activation function σ(x), it can be represented in the form of an ordinary differential equation (ODE) as follows: dσ dx(x) = σ(x)(1 −σ(x)).(3)
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 11 This formulation enables the use of ODE solvers within an affine arithmetic framework to compute guaranteed enclosures of the solutions while preserving the dependencies between individual neurons in the network. Supporting a broad spectrum of neural network architectures and expressive temporal logic specifications, the framework enables formal verification of practical NNCS scenarios. Comparative analysis with a Monte Carlo-based method highlights its precision and formal soundness. Figure 1: The 20-second simulation illustrates branching in the neural network output, with reachable tubes depicted in red. The axes indicate the position of the NN-controlled robot in meters. Branching arises from uncertainty in output classification. Left deviations around obstacles result in longer trajectories to the goal. Acknowledgement The authors acknowledge support from the French Interdisciplinary Center for Defense and Security (CIEDS) with the STARTS project.
Formal Verification of State and Temporal Properties of Neural Network-Controlled Systems Antoine Besset, Joris Tillet, and Julien Alexandre dit Sandretto Summary : Ensuring the safety of Neural Network Controlled Systems (NNCS) remains a major challenge due to the opacity of neural networks, especially when temporal properties are involved. By combining interval analysis with reachability techniques, we ensure compliance with spatial and temporal specifications formally expressed using Signal Temporal Logic (STL). The proposed framework provides a rigorous method for guaranteeing the correctness of uncertain dynamical systems controlled by a neural network. Cyber-physical systems Consider a continuous dynamical system modeled by the following differential equation: ˙ y(t) = ƒ(y(t), (t)), y(t)∈Rn, (t)∈W,(1) where y(t) is the state of the system, (t) is a bounded external input, and W⊆Rp is a compact set. For any initial state y0∈Rn and any measurable input :R+→W, the system admits a unique trajectory denoted by ξ(·, y0, ). In the presence of bounded uncertainty, the objective is to determine the set of possible trajectories over the interval [t0, T] . The set of reachable states at time t∈R+from an initial set Y0⊆Rnis defined as: Reacht(Y0,W) = {ξ(t, y0, )|y0∈Y0, (s)∈W,∀s∈[0, t]}.(2) Acontinuous-time representation on [tj, tj+1],SN−1 j=0[tj, tj+1] = [t0, T] , called a tube and denoted [˜ y](t) for t∈[t0, T] , is essential for preserving the set of all possible system behaviors [6]. This tube enables the analysis of the satisfaction of a temporal logic formula. Combining STL and reachability analysis To formally analyze the system behavior, we use interval analysis [5] and introduce a set - valued extension of predicates: ([ ˜ y], t)μ:= 1,if [˜ y](t)⊂Xμ, 0,if [˜ y](t)∩Xμ=∅, [0,1],otherwise. (3) Propagation in temporal logic is handled using Boolean intervals [2,7], e.g [3].: 0∧[0,1] = 0,0∨[0,1] = [0,1],1∧[0,1] = [0,1],1∨[0,1] = 1. Signal temporal logic We use the formalism of Signal Temporal Logic (STL) [4]: φ:=> | μ|¬φ|φ1∧φ2|φ1U[t1,t2]φ2 where the temporal operator U (Until) specifies that a property must hold until another becomes true within a given time interval. The operator F (Finally) expresses that a goal must be reached within a time window, while G (Globally) states that a property must hold throughout a time interval. Example: an automaton executing a periodic task. φ=G[0,5](⇒F[3.5,4.5])∧G[0,8](¬p)∧F[8,9]q •G[0,5](⇒F[3.5,4.5]) : Always on [0,5] s, if is reached then must follow within 3.5–4.5s. •G[0,8](¬p): The obstacle pmust never be encountered during the first 8s. •F[8,9]q: The stand-by zone qmust be reached between 8s and 9s. Tube ([ ˜ y], t)in blue and zones (Xμ) in red. Application : A robot controlled by a neural network A robot is controlled by a neural network, whose internal behavior is difficult to interpret. The presented methods provide formal guarantees that it reaches its goal, avoids obstacles, and does so within a given time bound, even in the presence of uncertainties [1]. The system employs a neural network to choose, from a set of motion primitives, the action that drives it toward the goal while avoiding obstacles. An example of specification could be: φ=¬CU[t1,t2]T. This means that no collision ( ¬C ) must occur until the target ( T ) is reached, within the given time horizon [t1, t2]. The reachable tube by the robot is shown in red, obstacles are in green and yellow, the target point is in purple. Set propagation in neural network To conduct reachability analysis of NNCS, activation functions such as the sigmoid can be expressed as ODEs, e.g. dσ d() = σ()(1−σ()). This allows ODE solvers to be combined with affine arithmetic, where uncertain quantities are represented as =0+1ϵ1+···+nϵn, ϵ∈[−1,1]. Using shared noise symbols preserves dependencies between neurons, enabling accurate error tracking and avoiding the overestimation of interval arithmetic. References [1] Antoine Besset, Julien Alexandre dit Sandretto, and Joris Tillet. Real-time guaranteed monitoring for a drone using interval analysis and signal temporal logic. In Proceedings of the 2025 IEEE/RSJ IROS, Hangzhou, China, 2025. IEEE. [2] Antoine Besset, Joris Tillet, and Julien Alexandre dit Sandretto. Uncertainty removal in verification of nonlinear systems against signal temporal logic via incremental reachability analysis. In Proceedings of the 64th IEEE CDC. IEEE, 2025. [3] Luc Jaulin, Michel Kieffer, Olivier Didrit, and Eric Walter. Applied interval analysis. In Luc Jaulin, Michel Kieffer, Olivier Didrit, and Eric Walter, editors, Applied Interval Analysis, pages 11–43. Springer, 2001. [4] Oded Maler and Dejan Nickovic. Monitoring temporal properties of continuous signals. In Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems, volume 3253, pages 152–166. Springer, 2004. [5] Ramon E. Moore. Interval Analysis. Series in Automatic Computation. Prentice Hall, 1966. [6] Julien Alexandre Dit Sandretto and Alexandre Chapoutot. Validated explicit and implicit runge-kutta methods. Reliable Computing, 22, 2016. Special issue devoted to material presented at SWIM 2015. [7] Joris Tillet, Antoine Besset, and Julien Alexandre Dit Sandretto. Guaranteed satisfaction of a signal temporal logic formula on tubes. Acta Cybernetica, 2025. Accepted, to appear. U2IS, ENSTA – Institut Polytechnique de Paris, 828 boulevard des Maréchaux, 91120 Palaiseau, France Projet STARTS Sécurité et Probabilité de Succès (SafeTy And pRobAbiliTy Success) starts.ensta-paris.fr
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 13 Exploiting the Impossible: Towards Resilience of Decision Making Against Misperceptions Martin Fr¨anzle, Paul Kr¨oger, and Anna Nienaber Carl von Ossietzky Universit¨at Oldenburg Foundations and Applications of Systems of Cyber-Physical Systems D-26111 Oldenburg, Germany {martin.fraenzle,paul.kroeger,anna.nienaber}@uni-oldenburg.de Keywords: Highly automated vehicles, learning-enabled cyber-physical systems, perception chain, decision making, robustification, symbolic-numeric computation Introduction One of the key challenges for safety-critical cyber-physical systems (CPSes) such as (highly) autonomous vehicles is decision making under the inevitable presence of uncertainties in environmental perception. Wrong control decisions within such systems may incur a substantial risk to life, health, or property. Achieving high confidence for guard conditions enabling safety-critical actions is thus crucial, even if they rely on uncertain percepts. However, the safety targets for, e.g., safety-critical manoeuvres of vehicles are typically orders of magnitude higher than, e.g., the statistical figures for the reliability of at least current learning-enabled object classification algorithms. The perception and decision chain consequently needs to incorporate mechanisms for substantially improving the confidence in critical guard conditions, i.e., to reduce the risk of erroneously performing an unsafe manoeuvre to a frequency considerably below the risk of individual misperceptions, e.g., misclassifications. These mechanisms should, however, not impede performance, i.e., they should retain liveness of the system in that they reduce the rate of erroneously admitting a safety-critical action drastically, yet do not significantly reduce the overall likelihood of permitting the respective action. We present a symbolic-numeric method that systematically rewrites critical guard conditions s.t. the resulting conditions are more resilient against misperceptions than the original conditions in that they are compatible with a given safety target, e.g. a societally accepted upper bound of the risk of erratically activating a safety-critical manoeuvre due to a false positive in a guard evaluation, while liveness in terms of the true positive rate of guard evaluation simultaneously is maximised. Approach 1 2 3 Figure 1: A traffic situation. Figure 1 illustrates a traffic situation in which the blue ego car (no. 1) shall overtake the orange car (no. 2) iff, first, the orange car is detected as an obstacle and, second, an overtaking manoeuvre is safe, i.e., iff there is no oncoming traffic (such
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 14 as the red car no. 3) with which car 1 could collide. In this –for the sake of conciseness over-simplified– example, the physical environment is partitioned into grid elements and the overtaking manoeuvre could be guarded by a complex Boolean condition describing conditions on the occupancy of the grid elements, which in turn is evaluated separately by a suitable but inherently uncertain object detection and classification algorithm, usually of machine-learning type. Such safe-guarding conditions are prone to induce unnecessary risk by demanding potentially unsafe overtaking manoeuvres as soon as a single grid element is (mis-) perceived. Incorporating environmental invariants can alleviate the problems described above: A car can usually neither occupy a single grid element only nor be distributed over non-adjacent grid elements. Our method exploits (formalisations of) such invariants. Given an invariant, we systematically rewrite a given guard by treating percepts not satisfying the invariant as don’t cares and generally remapping percepts s.t. the robustness of the guard condition against misperceptions increases. Robustness increases in that the rate of erratically evaluating the guard to be satisfied, i.e., its false-positive rate, is reduced to below a pre-defined threshold, while the true positive rate gets maximised in order to guarantee performance. We implemented and evaluated our approach by an algorithm that is based on reduced ordered binary decision diagrams (RoBDD) and akin to RoBDD don’t care optimisation. The algorithm is symbolic-numeric in that it attaches probabilities to the elements of an RoBDD and adjusts those numerically during the RoBDD operations underlying its don’t care optimisation. Acknowledgement This work was partially supported by the German Research Foundation (DFG) as part of PreCePT (FR 2715/6-1). References [1] M. Fr¨anzle and A. Hein. Safer Than Perception: Increasing Resilience of Automated Vehicles Against Misperception. Bridging the Gap Between AI and Reality (AISoLA 2023), LNCS 14129, 415–433, 2024. [2] Anna Nienaber. Increasing Resilience of Automated Vehicles Against Misperception. Bachelor’s thesis, Carl von Ossietzky Universit¨at Oldenburg, 2025. Unpublished.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 15 Computer-assisted proof of the simplicity of the second Dirichlet eigenvalue for non-equilateral triangles Ryoki Endo1, Xuefeng Liu2 1Graduate School of Science and Technology, Niigata University 8050 Ikarashi 2-no-cho, Nishi-ku, Niigata City, Niigata 950-2181, Japan [email protected] 2Department of Information and Sciences, Tokyo Woman’s Christian University 2-6-1 Zempukuji, Suginami-ku, Tokyo 167-8585, Japan [email protected] Keywords: Simplicity of eigenvalues, Dirichlet eigenvalue, Verified computation, Computer-assisted proof The rich relationship between Laplacian eigenvalues and shapes gave birth to the field of spectral geometry, which continues to attract researchers from various disciplines. In this talk, we provide a computer-assisted proof for a conjecture about Dirichlet eigenvalues posed by R. Laugesen and B. Siudeja in Henrot’s book “Shape Optimization and Spectral Theory” [1]: Conjecture 1 (Conjecture 6.47 of [1]).The second Dirichlet eigenvalue is simple on every non-equilateral triangle. The proof of this conjecture is given by two parts. In Part 1, we provided a partial result confirming the conjecture for the case of nearly degenerate triangles [2]: Theorem 1. The second Dirichlet eigenvalue is simple for every non-equilateral triangle with its minimum normalized height 1less than or equal to tan(π/60)/2. To achieve this, we derived explicit estimates for the k-th Dirichlet eigenvalues on the collapsing triangle: For s∈(−1,1) and t > 0, let T(s, t) be the triangular domain with vertices (−1,0),(1,0) and (s, t); see Figure 1. (s, t) −11 T(s, t) x Figure 1: Shape of triangle T(s, t) 1The minimum normalized height of a triangle is the height measured relative to its longest side, with the triangle scaled such that the longest side has unit length.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 16 Letting t0= tan(π/60)/2, ¯µk(s) 1 + t 2 3 0/(3π2)¯µk(s)≤t4 3λk(s, t)−π2 t2≤ˆµt0 k(s) (∀t∈(0, t0], k = 1,2,···),(1) where ˆµt0 k(s) is the k-th eigenvalue of a Schr¨odinger operator over a bounded interval, and ¯µk(s) is the k-th eigenvalue a Schr¨odinger operator on R. It is worth pointing out that the values or bounds of the involved eigenvalues are all computable by utilizing the recently developed methods for rigorous eigenvalue estimation [4]. The estimation for eigenvalues in (1) allows us to separate λ2(s, t) and λ3(s, t) for t∈(0, t0], confirming the simplicity of the second eigenvalue for nearly degenerate triangles. Part 2 completes the proof by covering the case of non-degenerate triangles [3]: Theorem 2. The second Dirichlet eigenvalue is simple for every non-equilateral triangle with its minimum normalized height greater than or equal to tan(π/60)/2. The methodology developed for this part provides a new way to stably compute eigenfunctions for clustered eigenvalues [5]. Acknowledgement Both authors are supported by Japan Society for the Promotion of Science. The first author is supported by JSPS KAKENHI Grant Number JP24KJ1170. The last author is supported by JSPS KAKENHI Grant Numbers JP22H00512, JP24K00538, JP21H00998 and JPJSBP120237407. References [1] Henrot, A. Shape optimization and spectral theory. De Gruyter Open Poland (2017). [2] Endo, R., Liu, X. The second Dirichlet eigenvalue is simple on every nonequilateral triangle Part I: Nearly degenerate triangles. Journal of Differential Equations 447, 113629 (2025). [3] Endo, R., Liu, X. The Second Dirichlet Eigenvalue is Simple on Every Non-equilateral Triangle, Part II: Nearly Equilateral Triangle. arXiv:2305.14063 (2024). [4] Liu, X. Guaranteed Computational Methods for Self-Adjoint Differential Eigenvalue Problems. Springer Nature (2024). [5] Endo, R., Liu, X. Stable computation of Laplacian eigenfunctions corresponding to clustered eigenvalues. arXiv:2506.07340 (2025).
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 17 Tuesday, September 23, 2025 A01-0-006 11:00–12:30 Regular Session B: Tools and Implementations 11:00–11:30 Ryoga Iwanami: Verry: an open-source package for verified computation written in Python 3 11:30–12:00 Pierre Filiol, Luc Jaulin, Theotime Bollengier and Jean-Christophe Le Lann: Hardware accelerated interval arithmetic for mobile robotics using RISC-V ISA extension 12:00–12:30 Jiˇr´ı Khun and Jan Schmidt: Towards Interval Arithmetic in TensorFlow: A Comparison of Approaches
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 18 Verry: an open-source package for verified computation written in Python 3 Ryoga Iwanami Waseda University Graduate School of Fundamental Science and Engineearing 3-4-1 Okubo, Shinjuku-ku, Tokyo, 169-8555, Japan [email protected] Keywords: delay differential equation, ordinary differential equation, Python In this talk, we present an overview of Verry [1], a verified computation library written in Python 3. Verry aims to provide a comprehensive implementation of various rigorous numerical algorithms for ODEs (e.g., [2], [3], [5], and [6]) and DDEs. We have released solvers based on [2] and [3] at this time. Since they solve the same problem, these algorithms have many common parts. We removed duplicate code by separating the algorithms into several routines. This separation also enables us to combine each routine depending on the problem. Here is an overview of the separation. Assume that a symbol put brackets around, like [x], always denotes an interval. Consider the initial value problem dy/dt=f(t, y) if t∈(t0, tbound), y∈[y0] if t=t0, where [y0]⊆RNis a non-empty bounded interval, and f(t, y) is a smooth function. A number of rigorous numerical algorithms for ODEs calculate the enclosure of the solution y= Φ(t, t0, y0) by the following iteration: 0. Initialize k= 0. 1. Set a coarse enclosure [pc k(h)] and verify that Φ(tk+h, tk, yk)∈[pc k(h)] holds for all h∈[0, tk+1 −tk] and yk∈[yk], where tk+1 is predefined or adaptively determined. Then refine [pc k(h)] into the tight enclosure [pk(h)]. 2. Calcurate [yk+1] such that Φ(tk+1, t0, y0)∈[yk+1] holds for all y0∈[y0]. 3. Increment k; then go to step 1 unless tkhas reached tbound. Note that one may obtain [yk+1] by evaluating an expression [pk(tk+1−tk)] directly; however, due to the wrapping effect [4], it induces an explosion of diam[yk]. We implemented separately steps 1 and 2 as abstract classes and showed the conditions that these subclasses must satisfy to cooperate with ODE solvers. For example, step 2 is implemented as an abstract class Tracker. The requirements that any xbeing an instance of Tracker must satisfy are as follows: 1. xcorresponds to some pair (c, S), where S⊆RNis star-shaped at c.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 19 2. x.sample() returns c, and x.hull() returns an interval vector containing S. 3. Given [A]⊆RN×Nand [b]⊆RN,x.update(A, b) updates (c, S) to (c′, S′) such that A(y−c) + b∈S′holds for all y∈S,A∈[A], and b∈[b]. Note that F(S)⊆S′holds if F∈C1(¯ S, RN), F(c)∈[b], and {DF(y)|y∈S} ⊆ [A]. These properties enable Tracker to track the trajectory of a given discrete dynamical system yn+1 =Fn(yn) (n= 0,1,2, . . . ). Hence, we can compute step 2 by applying Tracker to the system yn+1 = Φ(tn+1, tn, yn). Verry can also solve boundary value problems via the shooting method, and delay differential equations via the method of steps. We will show some numerical experiments. References [1] R. Iwanami,Verry, https://python-verry.github.io/verry (13 June 2025). [2] M. Kashiwagi,Power series arithmetic and its application to numerical validation, in Proceedings of the 1995 Symposium on Nonlinear Theory and its Applications, Las Vegas, NV, 1995, pp. 251–254. [3] R. J. Lohner,Enclosing the Solutions of Ordinary Initial and Boundary Value Problem, in Computerarithmetic, E. Kaucher, U. Kulisch, and Ch. Ullrich, eds., B. G. Teubner, Stuttgart, 1987, pp. 225–286. [4] R. J. Lohner,On the ubiquity of the wrapping effect in the computation of error bounds, in Perspectives on Enclosure Methods, U. W. Kulisch, R. J. Lohner, and A. Facius, eds., Springer-Verlag, Vienna, 2001, pp. 201–206. [5] N. S. Nedialkov and K. R. Jackson,An Interval Hermite–Obreschkoff Method for Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation, in Developments in Reliable Computing, T. Csendes, ed., Springer-Verlag, Dordrecht, 1999, pp. 289–310. [6] M. Neher, K. R. Jackson, and N. S. Nedialkov,On Taylor Model Based Integration of ODEs, SIAM J. Numer. Anal., 45 (2007), pp. 236–262, doi:10.1137/050638448.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 26 Interval and Set-Based Approaches for Control and State Estimation: Their Use for Offline and Online Purposes Andreas Rauh Carl von Ossietzky Universit¨at Oldenburg Distributed Control in Interconnected Systems D-26111 Oldenburg, Germany [email protected] Keywords: Interval methods, Set-based state estimation, Robust control, Optimization under uncertainty, Neural networks, Fuel cells Introduction Control and state estimation procedures need to be robust against imprecisely known parameters, uncertainty in initial conditions, and external disturbances. Interval methods and other set-based techniques form the basis for the implementation of powerful approaches that can be used to identify parameters of dynamic system models in the presence of the aforementioned uncertainties. Moreover, they are applicable to a verified feasibility and stability analysis of controllers and state estimators [1,3,7]. In addition to offline approaches for analysis, interval and set-based methods have also been developed in recent years which are allow to solve the associated design tasks and to implement reliable techniques that are applicable online. The latter approaches include online parameter adaptation techniques for nonlinear variablestructure controllers, interval observers, and fault diagnosis techniques [3,4,5,7]. In this talk, an overview of the methodological background will be presented, together with a review of practical applications for which interval and set-valued approaches have been employed successfully. Modeling, Parameter Identification, and Verified State Estimation Although, for example, many dynamic system models in (control) engineering, especially in the frame of thermo-fluidic applications, are described after a first-principle modeling by state equations that have certain monotonicity properties, other applications in the domain of mechanics as well as for electro-chemical energy storage may require specific changes of coordinates to obtain these properties. In the domains of parameter identification as well as state and disturbance estimation, the most important monotonicity property that allows for a simplification of the aforementioned tasks is the cooperativity of the state equations. As far as the application domains mentioned above are concerned, these properties originate from the conservation of mass or energy [4]. In such cases, a decoupling of lower and upper bounding systems — that enclose all possible state trajectories — can be obtained. This property does not only
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 27 allow for the simplification of the task of parameter identification but also for the implementation of real-time capable state estimation procedures. For systems with periodically recurring trajectories (and also disturbance profiles), recent investigations have shown that the corresponding procedures can also be extended to a learning-type technique. This technique especially allows for enhancing the bounds of estimated state trajectories in each successive execution of the same task and exploits a formulation that uses the iteration counter as a second independent dimension in addition to time [2]. Verified Control Implementation and Robust Model-Predictive Control On the basis of the set-based state estimates described in the previous section, realtime capable robust control implementations can be derived that prevent the violation of state constraints with certainty. Moreover, it is possible to implement robust predictive control laws in a similar manner. For the case of a nonlinear state feedback, interval extensions of sliding mode and backstepping control approaches have been published which allow for a guaranteed stabilization of the system dynamics and for a guaranteed prevention of overshooting certain thresholds for the state variables under constraints on the inputs and their respective variation rates. A practical application of this technique is the temperature control of a solid oxide fuel cell stack [3,7]. For the second class of controllers, a novel combination of set-based and neural network modeling was recently developed and integrated into a sensitivity-based predictive control scheme that maximizes the degree of fuel utilization of a fuel cell. The approach can be implemented for time-varying desired electric power profiles so that operating points stay within the region of Ohmic polarization, which is crucial for preventing accelerated aging of the fuel cell stack [5]. Combination of Set-Based and Stochastic Uncertainty Representations In the final part of this talk, a combination of stochastic and set-based (in this case, ellipsoidal) uncertainty representations will be considered. This approach allows, on the one hand, for a rigorous quantification of predefined confidence levels in stochastic state estimation procedures. On the other hand, it allows for handling nonlinearities in such a way that the previously mentioned tolerance bounds are definitely not determined in an overly optimistic manner [6]. References [1] E. Auer, L. Senkel, S. Kiel and A Rauh, Control-Oriented Models for SO Fuel Cells from the Angle of V&V: Analysis, Simplification Possibilities, Performance, Algorithms, 10, 140, 2017. [2] T. Chevet, A. Rauh, T.N. Dinh, J. Marzat, T. Ra¨ıssi, Robust Interval Observer for Systems Described by the Fornasini-Marchesini Second Model, IEEE Control Systems Letters, 6, 1940–1945, 2021.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 28 [3] N. Cont, W. Frenkel, J. Kersten, A. Rauh and H. Aschemann, Interval-Based Modeling of High-Temperature Fuel Cells for a Real-Time Control Implementation Under State Constraints, IFAC-PapersOnLine, 53, 12542–12547, 2020. [4] S. Ifqir, A. Rauh, J. Kersten, D. Ichalal, N. Ait-Oufroukh and S. Mammar. Interval Observer-Based Controller Design for Systems with State Constraints: Application to Solid Oxide Fuel Cells Stacks, Proceedings of the Intl. Conf. on Methods and Models in Automation and Robotics, MMAR 2019, Miedzyzdroje, Poland. [5] A. Rauh and E. Auer, Comparison of Stochastic and Interval-Based Modeling Approaches for the Online Optimization of the Fuel Efficiency of SOFC Systems, Proceedings of the 9th Intl. Conf. on Systems and Control (ICSC), 2021, Caen, France. [6] A. Rauh, T. Chevet, T.N. Dinh, J. Marzat, T. Ra¨ıssi, Robust Iterative Learning Observers Based on a Combination of Stochastic Estimation Schemes and Ellipsoidal Calculus, Proceedings of the 25th Intl. Conf. on Information Fusion (FUSION), 2022, Link¨oping, Sweden. [7] A. Rauh, L. Senkel and H. Aschemann, Interval-Based Sliding Mode Control Design for Solid Oxide Fuel Cells With State and Actuator Constraints, IEEE Transactions on Industrial Electronics, 62(8), 5208–5217, 2015.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 29 Green-Representable Solutions: Reformulating Suband Super-solution Theory for Poisson’s Equation Kazuaki Tanaka1, Ryoga Iwanami1, Kaname Matsue2and Hiroyuki Ochiai2 1Waseda University, Tokyo, Japan [email protected], [email protected] 2Kyushu University, Fukuoka, Japan {kmatsue, ochiai}@imi.kyushu-u.ac.jp Keywords: Poisson’s equation, Green’s function, Fundamental solution, Suband Super-solutions, Solution enclosure, Corner singularities Introduction We consider the weak solution u∈H1 0(Ω) of the following boundary value problem (−∆u=fin Ω, u= 0 on ∂Ω,(1) where Ω ⊂RN(N > 1) is a bounded polytopic domain, and fis a given function that satisfies f∈L1(Ω) for N= 1 and f∈Lp(Ω) for some p > N/2 when N≥2. The motivation of this research is to find upper and lower solutions that enclose the exact solution of this problem. However, with traditional definitions of upper and lower solutions, smoothness of functions is implicitly required, which made it impossible to represent upper and lower solutions using piecewise linear functions. To overcome this challenge, it is necessary to relax the conditions that upper and lower solutions must satisfy. Green-Representable Solutions We introduce a new framework based on fundamental solutions. For an evaluation point sint ∈Ω, we construct test functions of the form ϕsint (x) := aintΓ(sint, x) + Hsint (x),(2) where Γ is the fundamental solution satisfying −∆Γ(s, x) = δ(x−s), with explicit forms depending on dimension. aint is a non-zero coefficient, and Hsint is harmonic in some domain containing Ω. Definition 1 (Local Green-representability).A solution u∈H1 0(Ω) of the problem (1) is said to be Green-representable with respect to ϕsint (or simply ϕsint -Greenrepresentable for short) if, for a fixed sint, there exists a test function ϕsint of the form above such that aintu(sint) = ⟨f, ϕsint ⟩+Z∂Ω ∂u ∂nϕsint dγ. (3)
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 30 Definition 2 (Global Green-representability).A solution uof the problem (1) is said to be Green-representable with respect to mapping Φ : sint 7→ ϕsint (or simply Φ-Green-representable for short) if upossesses the regularity u∈H1 0(Ω) ∩W1,q(Ω) for some q > N, and there exists a mapping Φthat assigns to each point sint ∈Ωa test function ϕsint such that uis ϕsint -Green-representable. The following is our main theorem and its corollary for two-dimensional polygonal domains. Theorem 1 (Main Theorem).Let Ω⊂RN(N > 1) be a bounded N-dimensional polytopic domain and let f∈Lp(Ω) with p > N/2. Then, a solution u∈H1 0(Ω) of the problem (1) is locally Green-representable with respect to any fixed evaluation point in Ω. Corollary 1. Let Ω⊂R2be a bounded polygonal domain (possibly non-convex) and let f∈Lp(Ω) with p > 1. Then a solution u∈H1 0(Ω) of the problem (1) is globally Green-representable with respect to any mapping Φconstructed using fundamental solutions. Generalized Suband Super-Solutions Building on this representability, we define generalized suband super-solutions: Definition 3 (Green-representable super-solution).A function u∈W1,q(Ω) (q > N) is a Green-representable super-solution if there exists a nonnegative constant cand a mapping Φ : sint 7→ ϕsint such that ⟨∇u, ∇ϕsint ⟩ ≥ ⟨f, ϕsint ⟩+cZ∂Ω ∂ϕsint ∂n dγ (4) and u−c≥0on ∂Ω.(5) A corresponding definition applies to sub-solutions by reversing the inequalities. This generalization allows piecewise linear functions to serve as suband supersolutions—a capability not available in the classical framework. Theorem 2 (Comparison).Let ube a Green-representable solution with respect to mapping Φ, and assume that for all sint, we have ∂ϕsint ∂n ≤0on ∂Ω. Let uand ube a Green-representable sub-solution and super-solution, respectively, with respect to the same mapping Φ. Then, u−γsint ≤u≤u+γsint in Ω,(6) where γsint := 1 aint R∂Ω ∂u ∂n ϕsint dγ. In the talk, we will present numerical examples in one and two dimensions. In particular, for two-dimensional cases, we demonstrate pointwise evaluations in nonconvex domains where the solution lacks sufficient regularity.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 31 Verified Computation of All Positive Solutions to a H´enon-Type Boundary Value Problem Taisei Asai1, Kazuaki Tanaka2, Satoshi Tanaka3and Shin’ichi Oishi4 1Faculty of Science and Engineering, Waseda University 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan [email protected] 2Global Center for Science and Engineering, Waseda University 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan [email protected] 3Mathematical Institute, Tohoku University Aoba 6-3, Aramaki, Aoba-ku, Sendai 980-8578, Japan [email protected] 4Faculty of Science and Engineering, Waseda University 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan [email protected] Keywords: Verified numerical computation, All-solution search, Bifurcation analysis, H´enon-type equation We consider the H´enon-type equation, which is a two-point boundary value problem: (−u′′ = (|x|l+λ)up, x ∈(−1,1), u(−1) = u(1) = 0,(1) where the parameters land λsatisfy l≥0 and λ≥0, and the exponent psatisfies p > 1. In this context, an “even solution” refers to a function uthat is both a solution of (1) and an even function, satisfying u(−x) = u(x) for all x∈(−1,1). The case λ= 0 corresponds to the one-dimensional H´enon equation −u′′ =|x|lup, and problem (1) was introduced in [1] as part of a study on the symmetry of its solutions. In recent years, the H´enon-type equation has attracted attention due to the possibility of possessing more multiple solutions than the original H´enon equation, and the bifurcation structure of such solutions has been the subject of active research. According to the study [2], for fixed p > 1, the uniqueness of positive even solutions holds on most of the first quadrant (l, λ)⊂R2, and only a very narrow region remains as a candidate for the existence of multiple positive even solutions. However, while sufficient conditions for multiple solutions have been studied, the overall bifurcation structure—including the precise branching points—has not been fully understood. In this study, we rigorously determined the number of positive solutions of (1), including both even and non-even ones. The all-solution search is a method for rigorously identifying all solutions of a differential equation within a given parameter range. However, a fundamental difficulty in implementing this method arises from the fact that the set of parameters used to determine solutions (e.g., initial values) is noncompact. For example, to find all solutions of problem (1), one would, in principle,
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 32 need to explore the infinite range −∞ < u′(−1) <∞, which is computationally infeasible. To overcome this difficulty, we first established a priori estimates of the solutions. Proposition 1 gives an upper bound on the maximum value ∥u∥∞of any positive solution, while Propositions 2 and 3 provide explicit upper and lower bounds on the initial condition u′(−1), thereby reducing the search domain to a compact set. Based on these theoretical results, we constructed the initial value domain depending on the type of solution. For even solutions, we fixed u′(0) = 0 and varied u(0) within the bound given by Proposition 1. For general positive solutions (including non-even ones), we fixed u(−1) = 0 and varied u′(−1) within the bounds specified in Propositions 2 and 3. Here again, the upper bound from Proposition 1 plays a key role in computation: if the value of the function exceeds this bound during the numerical process, the trajectory can immediately be ruled out as a solution to (1), saving unnecessary computations. Using these settings, we carried out an efficient and exhaustive numerical search with the kv library [3]. Here, B(x, y) := R1 0tx−1(1 −t)y−1dt denotes the beta function. Proposition 1. Let ube a positive solution of (1) with l≥0,λ≥0, and p > 1. Then ∥u∥∞≤2B(l+ 1, p + 2) + λ p+ 2−1 p−1 . Proposition 2. Let ube a positive solution of (1) with l≥0,λ≥0, and p > 1. Then u′(−1) ≥1 2B(l+ 1,2) + λ 2−1 p−1 =1 21 (l+ 1)(l+ 2) +λ 2−1 p−1 . Proposition 3. Let ube a positive solution of (1) with l≥0,λ≥0, and p > 1. Then u′(−1) ≤s2(1 + λ) p+ 1 ·2p+1 2B(l+ 1, p + 2) + λ p+ 2−p+1 2(p−1) . Acknowledgement This work was supported by JSPS KAKENHI Grant Numbers JP23K19016, JP25K17304. References [1] S. Tanaka: Morse index and symmetry–breaking for positive solutions of one– dimensional H´enon type equations, Journal of Differential Equations, 255:7 (2013), 1709–1733. [2] N. Shioji, S. Tanaka, K. Watanabe: Multiple existence of positive even solutions for a two point boundary value problem on some very narrow possible parameter set, Journal of Mathematical Analysis and Applications, 513:1 (2022), 126182. [3] M. Kashiwagi: kv library, (2025). http://verifiedby.me/kv/
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 33 Semigroup approach for validating solutions to semilinear parabolic PDEs Akitoshi Takayasu1and Jean-Philippe Lessard2 1University of Tsukuba Institute of Systems and Information Engineering 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan [email protected] 2McGill University Department of Mathematics and Statistics 805 Sherbrooke West, Montreal, QC, H3A 0B9, Canada [email protected] Keywords: Semilinear parabolic PDEs, Initial value problems, Spectral methods, Semigroup theory, Computer-assisted proofs Introduction and contributions Recent advances in computer-assisted proofs for dynamical systems have been driven primarily by progress in the study of the global dynamics of infinite-dimensional problems, such as the rigorous construction of invariant objects, forward integration of time-dependent partial differential equations (PDEs), and the validation of connections between equilibria. In particular, solving the initial value problem (IVP) for PDEs has emerged as a central topic in this field. Over the past decades, several standard methodologies for rigorously integrating IVPs have been established, including the fully spectral approach [1], the autonomous semigroup approach [2], and the non-autonomous semigroup approach [3]. In this talk, we present our recent advances in the non-autonomous semigroup approach for semilinear parabolic PDEs, focusing on the validation of long-time existence of solutions and their asymptotic behavior. More precisely, we consider a class of IVP of the general form (ut= (λ0+λ1∆ + λ2∆2)u+ ∆pN(u), t > 0, x ∈Ω, u(0, x) = u0(x), x ∈Ω, where p∈ {0,1},Nis a polynomial satisfying both N(0) = 0 and its Fr´echet derivative DN(0) = 0, u0(x) is a given initial data. The parameters λ0,λ1,λ2are chosen so that the PDE is parabolic. The assumption of the PDE being semilinear implies that the degree pof the Laplacian in front of the nonlinear term Nis less than the one of the differential linear operator λ0+λ1∆ + λ2∆2. This form of PDEs covers wide variety of PDEs, such as the nonlinear heat, the Swift–Hohenberg, Cahn–Hilliard, Ohta–Kawasaki, phase-field-crystal (PFC), and Navier–Stokes, etc. Our approach reformulates the IVP as a zero-finding problem in a Banach space of time-dependent Fourier coefficients. A Newton-like operator is explicitly constructed
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 34 using the linearized evolution operator derived from semigroup theory. The inverse of the Fr´echet derivative is realized by the variation-of-constants formula, for which we develop rigorous numerical bounds using a decomposition into finite and infinite Fourier modes. A key feature of our approach is the explicit control of the evolution operator, which allows us to validate contraction properties of the Newton-like operator and thus prove the existence and local uniqueness of solutions in a neighborhood of numerical approximations. Furthermore, we extend this approach to a multi-step framework, enabling rigorous forward integration over long time intervals. In other words, this development relates to the challenge of efficiently controlling the wrapping effect in infinitedimensional problems, which is an essential problem in the field of interval analysis. By using the semigroup property, we rigorously control the evolution operator over multiple time intervals and thereby validate the contraction property of the Newtonlike operator over long time. This leads to more efficient rigorous integration of IVPs. In the talk, we demonstrate the effectiveness of this improved approach using the Swift–Hohenberg equation and the Ohta–Kawasaki equation, the latter of which includes derivatives in the nonlinear term. Finally, we discuss the parallelizability of the method, which makes it computationally efficient for high-dimensional systems. The non-autonomous semigroupbased rigorous integrator thus provides a unified and feasible framework for studying long-time dynamics of evolutionary PDEs through validated numerics. Acknowledgement AT is partially supported by the Top Runners in Strategy of Transborder Advanced Researches (TRiSTAR) program conducted as the Strategic Professional Development Program for Young Researchers by the MEXT and JSPS KAKENHI Grant Numbers 25K00922, 24K00538, and 22K03411. References [1] M. Cadiot, J.-P. Lessard. Recent advances about the rigorous integration of parabolic PDEs via fully spectral Fourier-Chebyshev expansions. arXiv:2502.20644 [math.AP], 2025. [2] J. B. van den Berg, M. Breden, R. Sheombarsing. Validated integration of semilinear parabolic PDEs. Numerische Mathematik, 156:1219–1287, 2024. [3] G. W. Duchesne, J.-P. Lessard, A. Takayasu. A Rigorous Integrator and Global Existence for Higher-Dimensional Semilinear Parabolic PDEs via Semigroup Theory. Journal of Scientific Computing, 102:62, 2025.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 35 Wednesday, September 24, 2025 A01-0-006 11:00–12:30 Regular Session B: Uncertainty Quantification 11:00–11:30 Olga Kosheleva and Vladik Kreinovich: For statistical analysis of big data, interval uncertainty is needed 11:30–12:00 Jahangir Alam, Ismail Hossain, Tausif Hossain, Md Nuruzzaman Sojib, Olga Kosheleva and Vladik Kreinovich: How to compare situations in which we measure different quantities with different uncertainty 12:00–12:30 Ekaterina Auer and Wolfram Luther: Towards Fair and Explainable Medical Risk Prediction Software via Dempster-Shafer Theory
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 42 Thursday, September 25, 2025 Library Auditorium 9:00–10:30 Plenary Lecture Tino Teige: Bringing Formal Methods from Academia to Real-World Applications in Industry: My Personal Two-Decades-Journey 10:30–11:00 Coffee Break 11:00–12:30 Regular Session A: Estimation 11:00–11:30 Quentin Brateau, Fabrice Le Bars and Luc Jaulin: Separator for the remoteness constraint 11:30–12:00 Damien Esnault, Simon Rohou, Fabrice Le Bars and Luc Jaulin: Computing Interval Detection-Probability Grids via Monte-Carlo method for Underwater Robotics 12:00–12:30 Marit Lahme and Andreas Rauh: Set-Based Identification of Characteristic Curves and Its Challenges for Real-Life Applications 12:30–14:00 Lunch Break – Food truck 14:30–16:00 Regular Session A: Linear Systems and Linear Algebra 14:30–15:00 Haruto Kijima and Takeshi Ogita: Fast Implementation of Interval Matrix Multiplication Using Infimum-Supremum Representation With SIMD Operations 15:00–15:30 Takeshi Terao, Yoshitaka Watanabe and Katsuhisa Ozaki: Verification of Singular Values under Oblique Inner Product Space for Matrices 15:30–16:00 Katsuhisa Ozaki and Toru Koizumi: Tight Enclosure of Matrix Multiplication using Fused Multiply-Add 16:00–16:30 Coffee Break 16:30–17:30 Regular Session A: Neural Networks 16:30–17:00 Attila Sz´asz and Bal´azs B´anhelyi: Parameter Robustness of Neural Networks 17:00–17:30 Tibor Csendes: Interval Based Verification of Adversarial Example Free Zones for Neural Networks 18:00–19:00 Meeting of the program committee including new member
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 43 Bringing Formal Methods from Academia to Real-World Applications in Industry – My Personal Two-Decades-Journey – Tino Teige BTC Embedded Systems AG Department of Innovation & Technology D-26135 Oldenburg, Germany [email protected] Keywords: Formal Methods, Model Checking, Constraint Solving Abstract In this invited talk I will give insights into my personal experience in bringing formal methods from academia to real-world applications in industry. I will first elaborate on the development of the constraint solver iSAT during my academic years at the University of Oldenburg from 2005 to 2012. After that I will address the challenges and use cases from real-world applications I was faced with when I moved to the industrial test tool provider BTC Embedded Systems in 2012. One particular focus of this talk will be on how powerful methods and tools from academia like iSAT have been successfully transferred from their academic environment to industrial software development projects.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 44 Separator for the remoteness constraint Quentin Brateau1, Fabrice Le Bars1and Luc Jaulin1 1Lab-STICC - ENSTA 2 rue Fran¸cois Verny, 29200 Brest, France [email protected],fabrice.le [email protected],[email protected] Keywords: Interval methods, State Estimation, Remoteness Introduction Interval analysis is well suited for solving state estimation problems, which can be useful to deal with uncertainty in dynamical systems [4, 5]. The remoteness constraint [6] constitutes a fundamental geometric relationship in robotics, defining the distance measure between sensors with measurement directivity and obstacles in the environment. This constraint is particularly critical for acoustic sensors, where the directivity pattern can exhibit substantial angular coverage [1], making precise characterization of sensor-obstacle spatial relationships essential for reliable navigation and mapping applications. We present an implementation of the separator [2] for the remoteness constraint. The proposed separator enables the complete characterization of the set of compatible sensor positions relative to obstacles, given distance measurements and sensor directivity constraints. Unlike previous approaches that were limited to inclusion testing [6], our implementation provides an efficient and modern way to characterize the feasible sensor positions. Main results Figure 1 shows a paving of separators [3] on the remoteness constraint relative to an obstacle segment shown in red, and the remoteness cone defined by vectors u1= −0.5−1and u2=0.5−1. The pink area represent the set of positions for the sensor compatible with the measured distance d, the blue areas represent the set of positions not compatible with the measurement, and the yellow area is the unknown area. Figure 1a shows the separator for the case d= [4,5], and Figure 1b shows the separator for the case d= [5,+∞]. This separator can be used to localize the set of possible positions for an underwater robot in a known environment such as a pool or a harbor. Figure 1c shows the state estimation of a robot in a pool performing cycles and taking two measurements in the pool. The set of starting positions for the cycle is well enclosed in the pink area. Acknowledgement This work has been supported by the French Government Defense procurement and technology agency (AID-2022850).
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 45 (a) d= [4,5] (b) d= [5,+∞] (c) Set of possible position of a robot performing cycles in a pool. Figure 1: Paving of separators on the remoteness constraint References [1] Yang Cong, Changjun Gu, Tao Zhang, and Yajun Gao. Underwater robot sensing technology: A survey. Fundamental Research, 1(3):337–345, 2021. [2] Luc Jaulin and Benoˆıt Desrochers. Introduction to the Algebra of Separators with Application to Path Planning. Engineering Applications of Artificial Intelligence, 33:141–147, August 2014. [3] Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter, Luc Jaulin, Michel Kieffer, Olivier Didrit, and ´ Eric Walter. Interval analysis. Springer, 2001. [4] Marit Lahme and Andreas Rauh. Set-valued approach for the online identification of the open-circuit voltage of lithium-ion batteries. Acta Cybernetica, 26:855–869, 11 2024. [5] Ghalia Nassreddine, Fahed Abdallah, and Thierry Denoeux. State estimation using interval analysis and belief function theory: Application to dynamic vehicle localization. IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, 40(5):1205–1218, 2010. [6] E. Seignez, M. Kieffer, A. Lambert, E. Walter, and T. Maurin. Experimental vehicle localization by bounded-error state estimation using interval analysis. In 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 1084–1089, 2005.
Separator for the remoteness constraint Quentin Brateau, Fabrice Le Bars, Luc Jaulin {quentin.brateau, fabrice.le_bars, luc.jaulin} @ ensta.fr ENSTA, 2 rue François Verny, 29200 Brest, CNRS, UMR 6285, Lab-STICC, Pôle IA & Océans, ROBEX Introduction The remoteness constraint is fundamental in mobile robotics. It represents the constraint between the distance measured by a sensor within a measurement cone and an obstacle segment. m ab h1 h2h u1 u2 Figure 1 – Remoteness of a segment [a,b]relative to a point mand two unit vectors u1and u2 Measured distance The remoteness constraint ensures that the measured distance dis on of {||ma||,||mb||,||mh||,||mh1||,||mh2||,+∞}, with ||mh|| =det(ab,am) ||ab|| (1) ∀p∈ {a,b},||mp|| =q(mx−px)2+ (my−py)2(2) ∀i∈ {1,2},||mhi|| =det(ab,am) det(ui,ab)(3) (4) Case hu1,abi ≥ 0∧hu2,abi ≥ 0 det(ab,am)≥0 det(u2,am)≥0 det(u1,am)≥0 det(u1,bm)≤0 ||mh1|| ||ma|| ab u1 u2 Figure 2 – Remoteness when hu1,abi ≥ 0∧hu2,abi ≥ 0 Case hu1,abi<0∧hu2,abi ≥ 0 det(ab,am)≥0 det(u2,am)≥0 hab,ami ≥ 0 hba,bmi ≥ 0 det(u1,bm)≤0 ||mb|| ||mh|| ||ma|| ab u1 u2 Figure 3 – Remoteness when hu1,abi<0∧hu2,abi ≥ 0 Case hu1,abi<0∧hu2,abi<0 det(ab,am)≥0 det(u1,bm)≤0 det(u2,bm)≤0 det(u2,am)≥0 ||mb|| ||mh2|| ab u1 u2 Figure 4 – Remoteness when hu1,abi<0∧hu2,abi<0 Results (a) u1=0.25 −1,u2=0.75 −1(b) u1=−2 −1,u2=−1 −1 (c) u1=0.75 −1,u2=0.75 −1 Figure 5 – Paving of the separator for d= [4,5] Application The remoteness constraint can be applied to the localization of an underwater robot in a known pool, using an echosounder to measure the distance to the pool walls, and by performing only two measurements. Figure 6 – Localization of an underwater robot in a known environment using remoteness constraint SCAN 2025, Oldenburg, Germany September 22 to 26, 2025
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 47 Computing Interval Detection-Probability Grids via Monte-Carlo method for Underwater Robotics D. Esnault1,2, S. Rohou2, F. Le Bars2and L. Jaulin2 1DGA Naval Systems (French Defense Ministry) Naval mine warfare department 29200 Brest, France [email protected] 2Lab-STICC, ENSTA ROBEX team 29200 Brest, France {simon.rohou,fabrice.le bars,luc.jaulin}@ensta.fr Keywords: Detection-Probability Grid, Monte-Carlo Method, Probabilities, Interval Methods and Underwater Robotics. Introduction Autonomous Underwater Vehicles (AUVs) are emerging as valuable assets in various domains due to their ability to operate autonomously in complex and hazardous environments [1]. To address the numerous sources of uncertainty—ranging from environmental disturbances (e.g., ocean currents, obstacles) to internal system limitations (e.g., sensor noise, state estimation errors)—AUVs are typically modeled and controlled using probabilistic approaches [2]. During mine-clearance missions, AUVs use onboard sensors to detect threats located on the seafloor. The area covered by these sensors during a mission can be represented as a set. Whether a target is detected then depends on whether its position falls within this set. However, due to the stochastic nature of the AUV’s motion, the coverage area itself is random. As a result, detection becomes a probabilistic event: each point on the seafloor has an associated detection probability that depends on the mission plan and the AUV’s behavior. To ensure mission effectiveness—especially when a list of potential mine positions is known—it is crucial to design the mission such that the detection probability for each target is close to 1. Detection probability grid A practical way to visualize the likelihood of detecting points on the seafloor is through a detection probability grid. This grid is a 2D matrix, where each cell corresponds to a specific point on the seafloor, and the associated value represents the probability of detection at that point [3]. We propose to compute an empirical estimation of this grid using a Monte-Carlo based method [3]. A large number of AUV simulations (trajectories) are generated based on a stochastic model. For each simulated trajectory, the area covered by the AUV’s
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 48 sensor is computed [4]. A 2D grid is then overlaid on the seafloor, where each point represents a potential location for detection. For each point and each simulation, a binary realization value is assigned: 1 if the point lies within the covered area, and 0 otherwise. By averaging these values across all simulations, we obtain an empirical estimate of the detection probability for each grid point. By applying this procedure across all points in the grid, we obtain a full detection probability map. Contribution In practice, it is often infeasible to perfectly compute the covered area for each AUV trajectory due to the high computational cost. Instead, we typically compute an approximation of the covered area with a guaranteed precision. This approximation introduces ambiguity: rather than dividing the space into only two categories (covered and not covered), the seafloor is partitioned into three complementary regions: definitely covered, definitely not covered, and ambiguous, where it is uncertain whether a point has been covered or not. Instead of a binary outcome, we must now work within a three-valued logic [5]. To handle this new framework, the Monte-Carlo method must be adapted accordingly. Recent work has proposed the use of interval-valued realizations to model such uncertainty [6]. In this approach, each realization is no longer a number but an interval: [1] when the point is definitely covered, [0] when it is definitely not, and [0,1] when the status is uncertain. By applying this interval-based Monte-Carlo method to each point in the 2D grid, the result is an interval detection probability grid, where each grid cell is associated with an interval of probability reflecting the uncertainty introduced by both uncertain factors and numerical approximations. References [1] R.B. Wynn, T. Le Bas, and al. Autonomous Underwater Vehicles (AUVs): Their past, present and future contributions to the advancement of marine geoscience. Marine Geology, vol. 352, no. 1, pp. 451–468, 2014. [2] P. Corke. Robotics, Vision and Control. Springer Nature, 2023. [3] S. Thrun, W. Burgard, and D. Fox. Probabilistic robotics. MIT Press, 2010. [4] B. Desrochers, and L. Jaulin. Computing a Guaranteed Approximation of the Zone Explored by a Robot. IEEE Transactions on Automatic Control, vol. 62, no. 1, pp. 425–430, 2016. [5] S.C. Kleene. Introduction to Metamathematics. Literary Licensing, LLC, 2012. [6] D. Esnault, S. Rohou, F. Le Bars, and L. Jaulin. Bounding the success probability of naval mine-clearance missions conducted by AUVs. to appear in the Proceedings of ”OCEANS 2025 - Brest”, Brest, France.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 49 Set-Based Identification of Characteristic Curves and Its Challenges for Real-Life Applications Marit Lahme and Andreas Rauh Carl von Ossietzky Universit¨at Oldenburg Distributed Control in Interconnected Systems D-26111 Oldenburg, Germany {marit.lahme,andreas.rauh}@uni-oldenburg.de Keywords: Identification, Interval Observer, Contractor, Battery System Introduction The knowledge of the relationship between different physical quantities of a system is essential for various applications, such as monitoring its current state, predicting its future behavior, and designing appropriate controllers. These relationships can, for instance, be graphically represented by characteristic curves. For example, the characteristic curve of an ideal Ohmic resistor, showing the current-voltage relationship, is a linear function in R2. In certain applications, it is not feasible to directly obtain the characteristic curves of interest by measuring the corresponding physical quantities. One example is the open-circuit voltage (OCV) characteristic of a lithium-ion battery, which represents the dependence of the OCV on the state of charge (SOC), a quantity that cannot be measured directly. This characteristic is particularly valuable for modeling the dynamic behavior of the battery and for detecting aging or degradation effects [1]. Typically, obtaining the OCV characteristic involves measuring the OCV during a controlled charging or discharging cycle over several hours, while simultaneously estimating the SOC. Our goal is to identify this characteristic during system operation; in earlier work, we therefore proposed an online set-based identification scheme for this purpose [1, 4]. In this presentation, we outline the identification scheme with a primary focus on the challenges encountered when applying it to real-life systems with the identification of the OCV characteristic of a lithium-ion battery as an example. This work does not present a novel contribution, but rather provides a comprehensive overview of the identification scheme in a real-life context, highlighting the associated challenges and selected solution approaches. Main results The proposed identification scheme is a two-stage procedure. Since the SOC cannot be measured directly, it has to be estimated. For this purpose, an interval observer is employed, which provides estimated lower and upper bounds for the SOC that are guaranteed to enclose the true value. Based on this estimate together with structural information resulting from the application of Kirchhoff’s voltage law to an equivalent circuit model the OCV can be computed, resulting in an interval enclosure for the true value of the OCV. The estimates for the OCV and SOC can be represented by interval boxes, which are the Cartesian products of their interval elements. During
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 50 battery charging/discharching, different interval boxes are obtained spanning different sections of the OCV-SOC characteristic. Intersecting overlapping interval boxes reduces the estimation uncertainty and enables the reconstruction of the characteristic curve [1]. Applying this identification scheme to real-life systems presents certain challenges, particularly in the state estimation part [2, 3, 5]. These challenges arise mainly from the structure of the system. The dynamic behavior of a lithium-ion battery can be captured by a quasi-linear state-space representation where the dynamics matrices depend on the SOC as one of the state variables. Furthermore, the state variables are scaled differently, and the output equation is nonlinear. A classical Luenberger observer approach (a variant of a linear model-based state reconstruction scheme) does not yield sufficiently accurate or feasible estimation results in this case and is therefore unsuitable [3]. Instead, an interval observer which provides more design degrees of freedom is used. However, determining suitable observer gains that lead to feasible and accurate estimation results is challenging for the case of the lithium-ion batteries. We present several techniques for systematically parameterizing the observer gains to obtain estimated lower and upper bounds with a small interval width [5]. Additionally, uncertainties have to be taken into account during the observer design, which include measurement noise, process noise and parametric uncertainty. We present different techniques to address this during the observer parameterization [3]. References [1] M. Lahme and A. Rauh. Set-Valued Approach for the Online Identification of the Open-Circuit Voltage of Lithium-Ion Batteries. Acta Cybernetica (Special Issue of SWIM 2022), 26(4):855–869, 2024. DOI: 10.14232/actacyb.301185. [2] M. Lahme and A. Rauh. Online Identification of the Open-Circuit Voltage Characteristic of Lithium-Ion Batteries with a Contractor-Based Procedure. Proc. of the 27th International Conference on Methods and Models in Automation and Robotics (MMAR), IEEE 39–44, 2023. DOI: 10.1109/MMAR58394.2023.10242524. [3] M. Lahme, A. Rauh, and G. Defresne. Interval Observer Design for an Uncertain Time-Varying Quasi-Linear System Model of Lithium-Ion Batteries. Proc. of the 2024 European Control Conference (ECC), IEEE 3696-3702, 2024. DOI: 10.23919/ECC64448.2024.10591102. [4] M. Lahme and A. Rauh. Combination of Stochastic State Estimation with Online Identification of the Open-Circuit Voltage of Lithium-Ion Batteries. Proc. of the 1st IFAC Workshop on Control of Complex Systems (COSY), IFACPapersOnLine 55(40):97–102, 2022. DOI: 10.1016/j.ifacol.2023.01.055. [5] M. Lahme and A. Rauh. Systematic Approach for Parameterizing TNL Interval Observers. 15th Summer Workshop on Interval Methods (SWIM 2024), Presentation, 2024. DOI: 10.13140/RG.2.2.33444.28800.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 51 Fast Implementation of Interval Matrix Multiplication Using Infimum-Supremum Representation With SIMD Operations Haruto Kijima1, Takeshi Ogita1 1Waseda University Graduate of Fundamental Science and Engineering 3-4-1 Okubo, Shinjuku-ku, Tokyo, Japan [email protected] Keywords: Interval Arithmetic, Matrix Multiplication, SIMD Operations Introduction Let Fbe a set of floating-point numbers. Let IF be a set of intervals whose endpoints are floating numbers. We propose a method of interval matrix multiplication using infimum-supremum representation. It is difficult to achieve high performance for interval matrix multiplication using the infimum-supremum representation. This is because in modern CPUs, comparing signs and switching rounding mode operations prevent efficient execution of floating-point operations. Therefore, Kn¨uppel [1] proposed a method of calculating a product of a point matrix and an interval matrix, taking advantage of the fact that a product of a floating-point scalar aand an interval vector [v] = [v,v] (v,v∈Fn) can be calculated as if a > 0then w=fl▽(a·v); w=fl△(a·v); else w=fl▽(a·v); w=fl△(a·v); However, since this is a level-1 BLAS operation, the execution performance is bounded by memory bandwidth. To overcome this, a fast algorithm of interval matrix multiplication was proposed using the midpoint-radius representation exploiting level-3 BLAS operations [2]. In the algorithm, the number of floating-point operations in the product of two interval matrices is 8n3. However, in the worst case, the resultant interval width is overestimated by 1.5 times. In response to this, Nguyen and Revol [3] proposed an algorithm that reduces the overestimation to 1.18, where the number of floating-point operations is 14n3. In addition, Rump [4] improved the algorithm, reducing the number of floating-point operations to 10n3. Proposed Method When implementing floating-point matrix multiplication, blocking technique is used for efficient computations. Here we focus on mb×nbregister blocking for interval cases. Let [A, A]∈IFmb×l,[B,B]∈IFl×nb,[C, C]∈IFmb×nbbe given. Consider computation of [C,C]←[A, A]·[B,B]+[C,C]. If the matrix multiplication is performed using rank-1 updates [5], the operation becomes [C, C]←[A∗1,A∗1]·[B1∗,B1∗] + ···[A∗l,A∗l]·[Bl∗,Bl∗]+[C, C]
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 58 When parameter perturbations are also taken into account, the robust training objective in Equation 1 is extended as shown in Equation 2, where λdefines the radius of allowable parameter perturbations: θ= arg min θ Ex,y[ max ˇ θ∈Bp1(θ,λ), ˇx∈Bp2(x,ϵ) Lˇ θ(ˇx, y)] (2) The Adversarial Weight Perturbation (AWP) algorithm approximates the nested maximization problem using an iterative, gradient-based optimization procedure, which leaves the model vulnerable to both inputand parameter-based attacks. In this work, the Adversarial Parameter Propagation (APP) algorithm was proposed to mitigate the shortcomings of AWP. Our method combines gradient-based search with interval arithmetic to compute a more accurate approximation of the inner maximization problem, thereby improving robustness against input and parameter perturbations. Results Several neural networks were trained based on the CNN7 architecture for the CIFAR10 image classification task, employing different values of the regularization parameter λ. As presented in Table 1, the proposed models show better performance on both clean and adversarial inputs compared to the AWP algorithm. The parameter robustness of our models was evaluated using the Adversarial Parameter Attack [4]. For smaller attack radii, our models demonstrated approximately 15% greater resistance, while for larger radii, we observed up to a 30% improvement in robustness. ϵ λ Algorithm Accuracy Robust Accuracy 2 255 0.01 APP 63.8% 49.63% AWP 60.68% 44.16% 0.02 APP 60.77% 48.18% AWP 58.1% 43.70% Table 1: Clean and robust accuracy on the CIFAR-10 test set References [1] Madry, A., Makelov, A., Schmidt, L., Tsipras, D. & Vladu, A. Towards Deep Learning Models Resistant to Adversarial Attacks. (2019) [2] Gowal, S., Dvijotham, K., Stanforth, R., Bunel, R., Qin, C., Uesato, J., Arandjelovic, R., Mann, T. & Kohli, P. On the Effectiveness of Interval Bound Propagation for Training Verifiably Robust Models. (2019) [3] Wu, D., Shu-Xia & Wang, Y. Adversarial Weight Perturbation Helps Robust Generalization. (2020) [4] Yu, L., Wang, Y. & Gao, X. Adversarial Parameter Attack on Deep Neural Networks. (2022)
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 59 Interval Based Verification of Adversarial Example Free Zones for Neural Networks Tibor Csendes1 1University of Szeged and University of Pannonia H-6720 Szeged, Hungary [email protected] Keywords: Adversarial Example, Artificial Intelligence, Interval Arithmetic, Verification Introduction Recent machine learning models are sensitive to adversarial input perturbation. That is, an attacker may easily mislead an otherwise well-performing image classification system by altering some pixels. It is quite challenging to prove that a network will have correct output when changing slightly some regions of the images. This is why only a few works targeted this problem. Although there are an increasing number of studies in this field, really reliable robustness evaluation is still an open issue. We will present some theoretical results on the dependency problem of interval arithmetic critical in interval based verification. Main results We investigate the overestimation amounts we can face while evaluating trained, fully connected feed forward networks with the ReLU activation function. We study the situation when the weights of such a network are given as real numbers, we fix an input (e.g. a picture), and we test how large intervals around the input values can be verified to result in the same classification we obtained for the real case. In the next theoretical investigations we assume that interval arithmetic is calculated in the precise way, i.e. we exclude the effect of outwards rounding. Proposition 1. For a fully connected feed forward standard artificial neural network the overestimation size w(F(X))−w(f(X)) of the inclusion function can be zero only if at least one of the following conditions are fulfilled: •all input intervals are of zero width: w(xi) = b−a= 0, •For each output, the weights associated with all input variables xihave the same sign, either all nonnegative or all nonpositive, and •all the final evaluation functions calculating the outputs of the network have negative arguments. These conditions are sufficient one by one, and a proper combination of them is also necessary.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 60 Consider now the question which are the major factors for the overestimation sizes in the same setting. Proposition 2. For a fully connected feed forward standard artificial neural network of kinput intervals, mneurons in each of the even number of nhidden layers, and all weights wibounded by |wi| ≤ W, the amount of overestimation w(F(X)) −w(f(X)) of the inclusion function of an output is not more than 2n/2mn/2WnPk i=1 w(Xi). Corollary 1. A direct consequence of Proposition 2 is that we can have the same amount of overestimation due to the dependency problem with decreasing the number of hidden layers while increasing the number of neurons in a layer and vice versa. The main consequence of our theoretical study is that we can control the amount of overestimation caused by the dependency effect of interval arithmetic by forcing advantageous parameters such as low absolute bound of weights, or minimizing the number of hidden layers – while keeping the expected level of precision and recall. Acknowledgement This research was supported by the project Extending the activities of the HUMATHS-IN Hungarian Industrial and Innovation Mathematical Service Network EFOP3.6.2-16-2017-00015, 2018-1.3.1-VKE-2018-00033, and by the Subprogramme for Linguistic Identification of Fake News and Pseudo-scientific Views, part of the Science for the Hungarian Language National Programme of the Hungarian Academy of Sciences (MTA). References [1] T. Csendes. An interval method for bounding level sets of parameter estimation problems. Computing, 41:75–86, 1989. [2] T. Csendes. Interval Based Verification of Adversarial Example Free Zones for Neural Networks – Dependency Problem. Annales Mathematicae et Informaticae, 60:19–26, 2024. [3] T. Csendes, N. Balogh, B. B´anhelyi, D. Zombori, R. T´oth, and I. Megyeri. Adversarial Example Free Zones for Specific Inputs and Neural Networks, Proc. of the 2020 ICAI, Eger, Hungary, https://ceur-ws.org/Vol-2650/paper9.pdf [4] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I.J. Goodfellow, and R. Fergus. Intriguing properties of neural networks. Proc. of the 2014 Int. Conf. on Learning Representations, doi: 10.48550/arXiv.1312.6199 [5] V. Tjeng, K. Xiao, and R. Tedrake. Evaluating Robustness of Neural Networks with Mixed Integer Programming. Proc. of the 2019 International Conference on Learning Representations, doi: 10.48550/arXiv.1711.07356 [6] D. Zombori, B. B´anhelyi, T. Csendes, I. Megyeri, and M. Jelasity. Fooling a complete neural network verifier. Proc. of the 2021 Int. Conf. on Learning Representations, https://openreview.net/forum?id=4IwieFS44l
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 61 Thursday, September 25, 2025 A01-0-006 11:00–12:30 Regular Session B: Uncertainty Quantification 11:00–11:30 Olga Kosheleva and Vladik Kreinovich: Inconsistencies in Fuzzy Estimations: Kaucher Arithmetic Naturally Appears 11:30–12:00 Miroslav Svitek, Olga Kosheleva and Vladik Kreinovich: Shapley Value Under Interval Uncertainty Revisited: Why Seemingly Natural Axiomatic Approach Is Not Fully Adequate 12:00–12:30 Luc Jaulin: A new wrapper for a reliable resolution of underdetermined nonlinear equations 14:30–16:00 Regular Session B: Optimization 14:30–15:00 Milan Hlad´ık: Linear Programming Problems with Absolute Values and Interval Uncertainty 15:00–15:30 Cyril Koteck´y and Milan Hlad´ık: Basis stability in interval quadratic programming 15:30–16:00 Christophe Jermann, Nathalie Revol and Christine Solnon: B&P algorithms for continuous constraint problems: a survey of branching strategies 16:30–17:30 Regular Session B: Dynamic Systems 16:30–17:00 Ramiz Dilji, Bernd Tibken, Robert Dehnert, Youping Fan, Regina Deisling and Laura Ackerschott: Estimation of the Domain of Attraction for Nonlinear Systems using the Bihari Inequality 17:00–17:30 Andreas Rauh and Marit Lahme: Observer-Based Approaches for a Verified Simulation and Pseudo State Estimation of Fractional Dynamic Systems
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 62 Inconsistencies in Fuzzy Estimations: Kaucher Arithmetic Naturally Appears Olga Kosheleva1and Vladik Kreinovich2 1Department of Teacher Education, University of Texas at El Paso El Paso, TX 79968, USA, [email protected] 2Department of Computer Science, University of Texas at El Paso El Paso, TX 79968, USA, [email protected] Keywords: Fuzzy logic, Interval-values fuzzy logic, Intuitionistic Fuzzy Logic, Interval Uncertainty, Kaucher Arithmetic Interval-values and Intuitionistic fuzzy logics: a brief reminder When experts describes how they solve tasks – e.g., how expert drivers control their cars – they usually use imprecise (“fuzzy”) words from natural language such as small. To describe this knowledge in precise computer-understandable terms, Lotfi Zadeh proposed to ask the expert to assign, to each possible value xof the corresponding property, a degree m∈[0,1] to which xsatisfies the property – e.g., to which xis small. He called this technique fuzzy. Experts often use logical connectives ∗– e.g., “and” and “or” – in describing their decisions. For example, a condition for a certain action may be that a car in front is close and that it brakes a little bit. To estimate the degree cof a statement A∗B based on degrees aand bof its component statements Aand B, Zadeh proposed to take c=f(a, b), where f: [0,1] ×[0,1] →[0,1] is a (non-strictly) increasing continuous function for which for a, b ∈ {0,1}, the value f(a, b) is the usual truth value of the corresponding logical operation. This approach led to many successes, but its representation of expert knowledge was not always perfectly adequate. Two ideas were proposed to make it more adequate. The first idea was to take into account that, just like an expert cannot describe the exact value of control – e.g., he/she only says “a little bit” – this same expert cannot meaningfully describe his/her degree of belief by a single number. The expert’s opinion would be described more adequately if we allow the expert to use the interval [m, m] of possible values. This is known as interval-valued fuzzy logic. In line with general interval techniques, once we know the degrees [a]=[a, a] and [b] = [b, b] of statements Aand B, it is reasonable to estimate the degree of A∗B as f([a],[b]) def ={a∗b:a∈[a], b ∈[b]}.Since f(a, b) is increasing, this leads to f([a],[b]) = [f(a, b), f(a, b)]. The second idea was to take into account that when the expert is not 100% sure that xis small, this means that he/she also has arguments that xis not small. So, in addition to the degree mto which the property is true, it makes sense to also ask for the degree m−to which this property is false. This is known as intuitionistic fuzzy logic. We assume that there is no inconsistency, so m+m−≤1. Then, when we have pairs (a, a−) and (b, b−) corresponding to Aand B, it is reasonable to apply f
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 63 to aand b, and to apply a de Morgan-dual operation g(a, b)def = 1 −f(1 −a, 1−b) to the values a−and b−. These two ideas have a different meaning, but from the purely mathematical viewpoint, they are equivalent: we can map an interval [x, x] to a pair (x, 1−x) and, vice versa, a pair (a, a−) to an interval [a, 1−a−]; then, all operations remain the same. Accounting for possible inconsistencies naturally leads to Kaucher arithmetic In some cases, there is an inconsistency between arguments for and against the same statement. In such cases, it makes sense to consider pairs (a, a−) for which a+a−>1; see, e.g., [1]. If we apply the above-mentioned transformation (a, a−)7→ [a, 1−a−] to such pairs, we get an “interval” [a, a] for which a > a. In interval computations community, such intervals are known as improper, with special arithmetic – first proposed by Kaucher – extending interval arithmetic to such intervals. As shown in [2], Kaucher arithmetic operations can be interpreted as follows: while the usual interval operations describe the set of all possible values f(a, b) when a∈[a] and b∈[b], Kaucher operations describe the intersection of all possible ranges of f(S) over all connected sets S⊆[a]×[b] whose projections to aand b-axis are exactly [a] and [b]. Our main result is that for every non-strictly increasing function f(a, b), the results of applying fto thus extended intuitionistic pairs correspond exactly to Kaucher arithmetic. Other relations between fuzzy and Kaucher arithmetic If we know the sets [a] and [b] of all possible values of aand b, then the set of all potentially possible values of a+bis the interval sum [a]+[b], while the Kaucher sum [a]+K[b] is the set of all definitely possible values. These two intervals are the simplest case of a family of embedded intervals – which is exactly what a fuzzy number is. References [1] L. V. Arshinsky, On te third forms of conjunction and disjunction in logic with vector semantics, Ontology of Designing, 2025, Vol. 15, No. 2, pp. 262–269 (in Russian). [2] V. Kreinovich, V. M. Nesterov, and N. A. Zheludeva, Interval methods that are guaranteed to underestimate (and the resulting new justification of Kaucher arithmetic), Reliable Computing, 1996, Vol. 2, No. 2, pp. 119–124.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 64 Shapley Value Under Interval Uncertainty Revisited: Why Seemingly Natural Axiomatic Approach Is Not Fully Adequate Miroslav Svitek1, Olga Kosheleva2, and Vladik Kreinovich3 1Faculty of Transportation Sciences, Czech Technical University in Prague 110 00 Prague 1, Czech Republic, [email protected] 2Department of Teacher Education, University of Texas at El Paso El Paso, TX 79968, USA, [email protected] 3Department of Computer Science, University of Texas at El Paso El Paso, TX 79968, USA, [email protected] Keywords: Shapley value, Machine Learning, Interval Uncertainty Shapley value: a brief reminder Many successes are due to collaboration, be it in manufacturing or in research. How to fairly divide the dividends between all nparticipants? For example, when we evaluate individual researchers, how to fairly distribute the overall points-for-paper between paper co-authors? In this division, it is reasonable to take into account what would be the productivity v(S) if only participants from the set S⊆Ndef ={1, . . . , n} worked together. The answer to this question was produced by the future Nobelist Lloyd Shapley. He formulated natural conditions: additivity; symmetry; and that a person who does not contribute anything, i.e., for whom v(S∪ {i}) = v(S) for all S, should not get anything. He proved that there is only one distribution scheme that satisfies these conditions, in which Person igets the amount xi(v) = Pa(|S|)·(v(S∪{i})−v(S)), where the sum is taken overall all sets Sfor which i∈ S,|S|denoted the number of elements in a set S, and a(m)def =m!·(n−m)! n!. This expression for xi(v) is known as the Shapley value. Lately, Shapley value has also been actively used in machine learning, to decide which of nfeatures used to make a decision are most important. In this case, v(S) is the effectiveness that we get when we only use features from the set S. Need for interval uncertainty In practice, we rarely know the exact values v(S). Often, we only know an interval [v](S)=[v(S), v(S)] that contains v(S). The agreement about division is usually decided before the project starts, in which case even the future value v(N) is not known exactly. In this case, a reasonable idea is to come up with intervals [x]i([v]) = [xi([v]), xi([v])]. Then we can use Hurwicz approach and make a distribution xi(v) = α·xi(v) + (1 −α)·xi(v),
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 65 where αis determined from the condition that the sum of these values should be equal to the overall amount v(N) – the overall monetary amount or the overall number of points for this particular paper. Current interval method and its limitation The paper [1] considers similar conditions to Shapley’s and shows that under these conditions, we should take, as bounds on xi([v]), the Shapley values corresponding to the functions v(S) and v(S). In many cases, this approach leads to reasonable results, but in other cases, it does not. For example, for n= 2, if v(S) = 0 for all S,v(∅) = v(∅) = v({x2}) = 0, and v({x1}) = v({1,2}) = 1, then Person 2 gets nothing, although it is possible, e.g., that the actual values are v({1,2}) = 1 and v({1}) = v({2}) = 0, in which case, due to symmetry, Person 2 should get exactly the same amount as Person 1. Analysis of the problem and resulting solution The reason for the above problem is that while the condition that v(S∪{i}) = v(S) for all Sindeed means that idid not contribute anything, but, as the above example shows, a similar interval equality [v](S∪{i})=[v](S) for all Sdoes not necessarily imply that Person iwas not contributing. So, a natural idea is to take, as [x]i([v]), the set of all possible values xi(v) for all functions vfor which v(S)∈[v](S) for all S. To find these intervals, let us take into account that the Shapley value formula can be reformulated as xi(v) = X S:i∈S a(|S|+ 1) ·v(S)−X S:i∈S a(|S|)·v(S). Thus, by using usual interval computations, we get: xi([v]) = X S:i∈S a(|S|+ 1) ·v(S)−X S:i∈S a(|S|)·v(S); xi([v]) = X S:i∈S a(|S|+ 1) ·v(S)−X S:i∈S a(|S|)·v(S). References [1] K. Autchariyapanikul, O. Kosheleva, and V. Kreinovich, “Shapley value under interval uncertainty and partial information”, In: V. Kreinovich, W. Yamaka, and S. Leurcharusmee (eds.), Data Science for Econometrics and Related Topics, Springer, Cham, Switzerland, to appear.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 66 A new wrapper for a reliable resolution of underdetermined nonlinear equations Luc Jaulin ENSTA, Lab-STICC, France [email protected] Keywords: Nonlinear equations, Underdetermined system, Interval Methods Introduction This paper introduces a new wrapper called a buche, the French name for log (think of logs made from a straight trunk obliquely and bluntly cut with an axe). Buches are used to enclose a part of the solution set defined by nonlinear equations. We show that buches, combined with interval methods [2], may allow us to obtain a better accuracy for the approximation with less computations. Notion of buche The buche associated with a box [x]⊂Rn, a matrix A, a vector band the inflation rate ρis the set ⟨x⟩defined by ⟨x⟩=⟨[x],A,b, ρ⟩ ={x∈[x],∃p,Ap =band ∥x−p∥< ρ}.(1) An illustration is given by the figure below. The quantity ρ= rad(⟨x⟩) is called the radius of the buche ⟨x⟩. The affine space Ap =bis called a flat. Our motivation for using buches is to have the following properties •The box [x] in the structure of the buche will allow us to build a nonoverlapping covering of X. This is not the case for zonotopes [4], [1]. •A buche can easily be bisected, contrary to ellipsoids [3]. •The axis-aligned projection is easy with buches, contrary to polyhedrons. •A first order approximation is possible, contrary to boxes.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 67 The buche (green) of the picture corresponds to the intersection between the box [x] and a cylinder Contribution Buches will be used to represent the solution set of an underdetermined set of nonlinear equations. It will be shown that the use of buches makes it possible to increase the accuracy of the approximation of the solution set compared to classical interval techniques. References [1] C. Combastel. A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes. In 44th IEEE Conference on Decision and Control, 2005. [2] R. Moore. Methods and Applications of Interval Analysis. Society for Industrial and Applied Mathematics, jan 1979. [3] A. Rauh, J. Soueidan, S. Rohou, and L. Jaulin. Experimental validation of an ellipsoidal state estimation procedure for a magnetic levitation system. In IFAC World Congress, 2023, Yokohama, Japan. pp.8494-8499. [4] J. Wan. Computationally reliable approaches of contractive model predictive control for discrete-time systems. PhD dissertation, Universitat de Girona, Girona, Spain, 2007.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 74 Estimation of the Domain of Attraction for Nonlinear Systems using the Bihari Inequality Ramiz Dilji, Bernd Tibken, Robert Dehnert, Youping Fan, Regina Deisling and Laura Ackerschott University of Wuppertal Institute of Automatic Control D-42119 Wuppertal, Germany {dilji,tibken,dehnert,yofan,deisling,laura.ackerschott}@uni-wuppertal.de Keywords: stability analysis, nonlinear systems, integral inequalities Introduction In this contribution, we apply a long and well known integral inequality to the stability analysis of nonlinear systems. The main focus is on nonlinear systems with the state space representation ˙x=Ax +g(x),x(0) = x0,(1) where x∈Rnrepresents the state vector, A∈Rn×nis assumed to be Hurwitz, g represents the nonlinear part of the system starting with quadratic terms, and x0is the initial condition, respectively. Note that in this work, matrices appear as bold italic uppercase letters, column vectors as bold italic lowercase letters, and scalars as italic lowercase letters. We assume that the solutions of system (1) exist for all t≥0. Due to the eigenvalue condition of A, the origin is asympotically stable, and the problem of estimating the domain of attraction arises. Using the well known method of Lyapunov is the main approach to tackle this problem. In this abstract, we will use a different approach based on the following theorem of Bihari [1]. Theorem 1. Let ube a non-negative, continuous function satisfying u(t)≤α+Zt 0 f(s)w(u(s)) ds, t ∈[0,∞),(2) where α≥0is a constant, fis a non-negative continuous function and wis a continuous non-decreasing function with w(u)>0for all u > 0. Then it holds that u(t)≤G−1G(α) + Zt 0 f(s)ds, t ∈[0, T],(3) where G−1is the inverse function of G(x) = Zx x0 dy w(y), x ≥0, x0>0, and Tis determined such that the expression inside G−1in inequality (3) remains within the domain of G−1for all t∈[0, T].
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 75 Application to the Nonlinear State Space Representation The nonlinear system (1) is equivalent to the following integral equation x(t) = eAtx0+Zt 0 eA(t−τ)g(x(τ)) dτ. Taking the norm on both sides, we obtain ∥x(t)∥ ≤ M∥x(0)∥e−βt +MC e−βt Zt 0 e−βτ u(τ)2dτ, (4) where M > 0, β > 0 and C > 0 are constants. Here, Mand βarise from the bound ∥eAt∥ ≤ M e−βt, and Cprovides an upper bound for the nonlinear term according to ∥g(x(τ))∥ ≤ C∥x(τ)∥2. The constant βcorresponds to the absolute value of the real part of the rightmost eigenvalue of the matrix A. This choice of βensures exponential decay of the linear part. Ideally, βshould be as large as possible, while Mand Cshould be as small as possible to ensure tight estimates. Introducing the substitutions u(t) = eβt ∥x(t)∥, α =M∥x(0)∥, f (τ) = MC e−βτ ,and w(y) = y2, we can now rewrite system (1) into inequality (2). This allows the use of inequality (3), which provides u(t)≤1 1 M∥x(0)∥−MC β(1 −e−βt). For this inequality to be valid, the denominator on the right-hand side must remain positive. Taking the limit as t→ ∞, we obtain the condition ∥x(0)∥ ≤ β M2C, which provides an estimation of the domain of attraction for system (1). This estimate can, under certain circumstances, be larger than the one obtained using Lyapunov’s method, which corresponds to a more precise estimate of the domain of attraction. This is caused by the fact that any matrix norm can be used to build the inequality (4). By contrast, Lyapunov’s method usually requires a symmetric, positive-definite matrix to define the norm. This can result in conservative estimates of the domain of attraction, whereas the novel approach can provide a less conservative estimate. References [1] I. Bihari. A generalization of a lemma of Bellman and its application to uniqueness problems of differential equations. Acta Mathematica Academiae Scientiarum Hungaricae 7, 81–94 (1956). https://doi.org/10.1007/BF02022967
0 41,40 41,40 Hintergrundfarbe ändern: Platzhalter auswählendann über Menü: Format | Fülleffekt | Benutzerdefinierte Farbe auswählen 21,00 21,00 Hintergrundfarbe ändern: Platzhalter auswählendann über Menü: Format | Fülleffekt | Benutzerdefinierte Farbe auswählen 58,80 58,80 Folie in Ursprungsform bringen über Menu: Start | Folien | Zurücksetzen Weitere Formatierungen über Menu: Start | Absatz | Listenebene erhöhen/verringern Unbedingt angeben: Quellenangabe Foto (direkt in das Feld oder über Menü: Einfügen | Kopfund Fußzeile) Logo einfügen über den Bild Button (Grafik aus Datei einfügen) –ggf. Größe anpassen – wenn nicht benötigt, Platzhalter löschen School of Electrical, Information and Media Engineering Estimation of the Domain of Attraction for Nonlinear Systems using the Bihari Inequality Ramiz Dilji, Bernd Tibken, Robert Dehnert, Youping Fan, Regina Deisling, Laura Ackerschott 1Problem Statement •Stability analysis of nonlinear systems with state-space representation : state vector : system matrix, assumed to be Hurwitz : nonlinear part, starting with quadratic terms : initial condition •Due to eigenvalue condition of , the origin is asymptotically stable •Goal: Estimation of the domain of attraction for the origin •Classical approach: Lyapunov’s method →requires a symmetric, positive-definite matrix to define the norm →can lead to conservative estimates of the domain of attraction for nonlinear systems References Rainer-Gruenter-Str. 21 42119 Wuppertal, Germany University of Wuppertal Chair of Automatic Control Prof. Dr.-Ing. Bernd Tibken 2Theorem of Bihari Let be a non-negative, continuous function satisfying where is a constant, is a non-negative continuous function and is a continuous non-decreasing function with for all . Then it holds where is the inverse function of and is determined such that the expression inside xxx in inequality remains within the domain of xfor all xxxx . = 3Application •Equivalent integral equation of nonlinear system : •Estimation of the norm: : constants •Taking the norm: •Rewrite System into inequality by introducing the xlifollowing substitutions •The use of inequality results in •For this inequality to be valid, the denominator must llliremain positive. Taking the limit as leads to lllithe condition •Enables estimation of the domain of attraction for lxsystem 4Example System: 5Conclusion / Outlook •Novel approach •Application of Bihari’s Theorem for estimation of the domain of attraction •Reformulates the problem into an integral inequality •Provides an explicit bound on system trajectories •Any matrix norm that satisfies the triangle inequality is acceptable →a positive-definite matrix is not required →here, for example, the Euclidean matrix norm was used for the presented example •Can provide a less conservative estimate than Lyapunov’s method •Easier to compute, no hard optimization problem needs to be solved 6References [1] I. Bihari. A generalization of a lemma of Bellman and its application to uniqueness problems of differential equations. Acta Mathematica Academiae Scientiarum Hungaricae 7, 81-94 (1956). https://doi.org/10.1007/BF02022967
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 77 Observer-Based Approaches for a Verified Simulation and Pseudo State Estimation of Fractional Dynamic Systems Andreas Rauh and Marit Lahme Carl von Ossietzky Universit¨at Oldenburg Distributed Control in Interconnected Systems D-26111 Oldenburg, Germany {andreas.rauh,marit.lahme}@uni-oldenburg.de Keywords: Fractional dynamic systems, Finite memory approximation, Verified state estimation Introduction Fractional system models [1, 3] have gained in importance during recent years. They account for non-standard dynamics with long-term memory effects. Such kind of dynamics can be found, for example, in electrochemical energy converters such as batteries and fuel cells. Despite capturing infinite horizon memory properties, the numerical evaluation may be complicated by a continuous increase in memory and computing time if no appropriate countermeasures are taken. This is especially critical in the frame of (pseudo) state estimation, where a periodic integrator reset takes place to limit both memory demand and computing times in real-time applications. Additionally, such resets occur when predictor–corrector state estimators are employed for systems with bounded uncertainty. Solution Approach This contribution provides an extension of recent work [4], in which the authors have proposed a novel method that allows for estimating errors in a set-based form that result from truncating the memory of fractional system simulators to a finite length. The general idea is based on error bounds published in [3] for Riemann-Liouville fractional differential equations. It is then extended by an interval observer similar to the one in [1] after including the approximation errors by means of integrator disturbance models. A verified enclosure of the solution sets becomes possible after converting the point-valued observer of Theorem 1 into a set-based formulation, cf. [4]. Theorem 1 ([4] Point-Valued Observer for Truncation Errors).The fractional differential equation model (derivative order 0< ν ≤1, initialized at t0+T) "t0+TD(ν) tˆ z(t) t0+TD(ν) tˆ µ(t)#=f(ˆ z(t)) + ˆ µ(t) 0+H·(ym(t)−ˆ y(t)) (1) with ˆ z(t) = 0for t≤t0+Tis an observer for state (ˆ z) and additive truncation error reconstruction (ˆ µ) after an integrator reset at the point t0+Twith the sensor data ym(t)and the associated measurement model ˆ y(t) = h(ˆ z(t)). For stability, the gain Hneeds to be chosen so that the error dynamics are asymptotically stable.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 78 Fig. 1 gives an example for the observer-based state reconstruction and truncation error estimation for the system model x(0.5)(t) = −x(t)+u(t), z(t) = x(t)−x(t0+T), u(t) = 1 with the fractional differentiation order 0.5. It further contains a comparison between the integrator resetting with and without the proposed observer approach (obs.) and the true pseudo state evolution. (a) Simulation of x(t). (b) Illustration of the integrator resetting. (c) Approximation error of x(t). (d) Truncation error µ(t). Figure 1: Illustration of the observer-based state and truncation error estimation. References [1] E. Hildebrandt, J. Kersten, A. Rauh, and H. Aschemann. Robust Interval Observer Design for Fractional-Order Models with Applications to State Estimation of Batteries, Proc. of the 21st IFAC World Congress, Berlin, Germany, 2020. [2] G. Bel Haj Frej, R. Malti, M. Aoun, and T. Ra¨ıssi. Fractional Interval Observers and Initialization of Fractional Systems. Communications in Nonlinear Science and Numerical Simulation, 82:105030, 2020. [3] I. Podlubny. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Mathematics in Science and Engineering, Academic Press, London, 1999. [4] A. Rauh and M. Lahme. A Finite Memory Approach Applied to Verified Pseudo State Estimation of Fractional Models of Lithium-Ion Batteries. IFACPapersOnLine, 12th IFAC Conference on Fractional Differentiation and its Applications, 58:185–190, Bordeaux, France, 2024.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 79 Friday, September 26, 2025 Library Auditorium 9:00–10:30 Plenary Lecture Christoph Matheja: Automated Verification of Discrete Probabilistic Programs 10:30–11:00 Coffee Break 11:00–12:30 Regular Session A: Dynamic Systems 11:00–11:30 Robert Szczelina, Anna Gierzkiewicz and Jakub Kural: Investigating chaos in Delay Differential Equations with rigorous numerical methods 11:30–12:00 Jakub Kural, Anna Gierzkiewicz and Robert Szczelina: Computer assisted proof of existence of periodic solutions to ENSO delay differential equation model 12:00–12:30 Th´eo Le Terrier, Marie Babel and Vincent Drevelle: Ultra-wideband Based Smart Wheelchair Pose Estimation using Interval Analysis 12:30–14:00 Lunch Break – Food truck 14:00–15:30 Regular Session A: Dynamic Systems 14:00–14:30 Andreas Rauh and Friederike Bruns: Set-Based Contracts for Systematic Controller Tuning in Interconnected Dynamic Systems 14:30–15:00 Anna Gierzkiewicz, Maciej Capinski and Pau Martin: Oscillating orbits in the Sitnikov model: equal masses case 15:00–15:30 Mohamed Fnadi and R´egis Lherbier: Interval Particle Filter for LiDAR-Based Object Tracking 16:00–. . . Closing of SCAN 2025
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 80 Automated Verification of Discrete Probabilistic Programs Christoph Matheja Carl von Ossietzky Universit¨at Oldenburg Theory of Correct Systems D-26129 Oldenburg, Germany [email protected] Keywords: probabilistic programs, deductive verification, formal methods Introduction Probabilistic programs are ordinary computer programs with the ability to base decisions on samples drawn from probability distributions. They appear as implementations of randomized algorithms and communication protocols, as well as descriptions of physical and statistical models (cf. [1] for an overview). Common questions in the analysis of probabilistic programs concern quantifying their expected behavior, e.g. how large is the expected runtime of an algorithm, the expected number of retransmissions of a network protocol, or the probability that a particle reaches its destination? Writing correct probabilistic programs is notoriously hard, arguably even harder than ordinary software development [3]. Over the last 15+ years, verification techniques for probabilistic programs have thus received much attention. By now, there exists a plethora of proof techniques for quantifying, amongst others, the termination probability or expected runtimes of such programs. To enable reasoning about the correctness of feature-rich probabilistic programs, those techniques must be adapted and combined. However, many of the existing verification techniques are inspired by different fields in computer science, mathematics, and engineering, such as control theory, program logics, probabilistic model checking, probability theory, and domain theory. Comparing—let alone combining—those techniques can thus be non-trivial (cf. [8]). Modern program verifiers often have a front-end that translates a program and its specification into an intermediate language, such as Boogie [6]. Such intermediate languages enable the encoding of complex verification techniques, while allowing for the separate development of efficient back-ends, e.g. verification condition generators or symbolic execution engines. In the same spirit, HeyVL [7] is a recently developed quantitative intermediate verification language, which aims to enable researchers to (i) prototype and automate new verification techniques for probabilistic programs, (ii) combine those techniques, and (iii) benefit from improvements to common backends. It is part of the caesar automated verification infrastructure2, which has been applied successfully to analyze probabilistic programs with different techniques. Figures 1 and 2 depict examples of probabilistic programs together with a description of the verified property and the technique encoded in HeyVL. 2https://www.caesarverifier.org
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 81 while (1 < i){ n:= i; while (0 < n){ d:= flip(0.5); i:= i−d;n:= n−1 } } Figure 1: Model of Rabin’s Mutal Exclusion Protocol [4]. Property: the probability to select exactly one process (i.e. i= 1) plus the probability of nontermination is at least 2 3if n≥2 holds initially. Verified by encoding weakest liberal preexpectations [5]. while (0 < x){ i:= N+ 1; while (0 < x < i){ i:= unif(1, N) };x:= x−1 } Figure 2: Model of the Coupon Collector’s Problem. Property: the expected number of loop iterations is bounded from above by N·HN, where HNis the N-th harmonic number. Verified by encoding the expected runtime calculus [2]. Outline This lecture will provide an overview of verification techniques for discrete probabilistic programs. Along the way, we will demonstrate how those techniques can be automated using the HeyVL intermediate language and caesar. References [1] G. Barthe, J.-P. Katoen, and A. Silva. Foundations of Probabilistic Programming. Cambridge University Press, 2020. [2] B. L. Kaminski, J.-P. Katoen, C. Matheja, and F. Olmedo. Weakest precondition reasoning for expected runtimes of randomized algorithms. Journal of the ACM, 65(5), 1-68. 2018. [3] B. L. Kaminski, J.-P. Katoen, and C. Matheja On the Hardness of Analyzing Probabilistic Programs. Acta Informatica 56(3), 255–285, 2019. [4] E. Kushilevitz and M. O. Rabin. Randomized Mutual Exclusion Algorithms Revisited. PODC. 1992. [5] A. McIver and C. Morgan. Abstraction, Refinement and Proof for Probabilistic Systems. Monographs in Computer Science., Springer, 2005. [6] K. R. M. Leino. This is Boogie 2. 2008. [7] P. Schr¨oer, K. Batz, B. L. Kaminski, J.-P. Katoen, and C. Matheja. A Deductive Verification Infrastructure for Probabilistic Programs. Proc. ACM Program. Lang. 7(OOPSLA2): 2052–2082, 2023. [8] T. Takisaka, Y. Oyabu, and N. Urabe, and I. Hasuo. Ranking and Repulsing Supermartingales for Reachability in Randomized Programs. ACM Trans. Program. Lang. Syst. 43(2), 5:1–5:46 2021.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 82 Investigating chaos in Delay Differential Equations with rigorous numerical methods Anna Gierzkiewicz1, Jakub Kural1, Robert Szczelina1,2 1Jagiellonian University Institute of Computer Science and Computational Mathematics Go l¸ebia 24 street, 31-007 Krak´ow, Poland 2{robert.szczelina}@uj.edu.pl Keywords: Delay Differential Equations, Pseudospectral method, Interval Newton Operator, Infinite-dimensional Dynamical Systems Abstract We are studying the Cauchy problem for Delay Differential Equations with constant delays of the following form: x′(t) = f(x(t), x(t−τ)) t≥0 x(t) = ψ(t)t∈[−τ, 0], with τ > 0 fixed, x∈Rd,ψ∈ C0([−τ, 0],Rd) =: C. In recent years, a significant effort was made to study and prove various dynamical phenomena in such systems such as periodic motions, with a particular attention given to proving chaos in canonical examples such as the Mackey–Glass equation [9]. Several of the techniques require rigorous numerical computations of a considerable scale, for instance, see [1, 3, 4, 7, 8, 10, 11] and references therein. However, the infinitedimensional nature of the DDEs presents some challenges, such as the complicated setup, finding good approximations, and carrying out the computer assisted proofs. For several years we have been working on methods for proving complicated dynamics (for example, periodic orbits and symbolic dynamics in DDEs) using rigorous, forward in time integration of the DDE in the (subspace of) the phasespace C[5, 11]. However, our current method uses a very high dimensional projection of the solutions, and thus is unfeasible for numerical investigations, especially when preparing data for computer assisted proofs, for instance, when looking for initial sets in computer assisted proofs. To address this issue we propose another approach: by using a pseudospectral approximation [2] we reduce the DDE to a finite-dimensional system of Ordinary Differential Equations (ODEs) while preserving numerically observed dynamical features of the original system. Due to the low-dimensionality of the resulting approximation, the computations are less demanding and can be done using known tools, such as CAPD rigorous ODE solvers [1], to efficiently verify some dynamical phenomena that closely mirror those of the full system. We present some rigorous results for both DDE and the approximating projection to ODE and we discuss some problems that arise in this approximation. Based on the experiments, we believe those complications are intrinsic to the nature of DDEs and solving them in the reduced system might guide the search for the computer assisted proof of chaotic (symbolic) dynamics in the full DDE system.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 83 Acknowledgement The author would like to acknowledge the support of Polish National Science Center (NCN) grant no. 2023/49/B/ST6/02801 and Polish NAWA Bekker project no. BPN/BEK/2023/1/00170. References [1] F.A. Bartha, T. Krisztin, and A. V´ıgh. Stable periodic orbits for the Mackey–Glass equation. J Diff.l Eq., Vol. 296 pp. 15–49, 2021. [2] D. Breda, O. Diekmann, M. Gyllenberg, F. Scarabel, and R. Vermiglio. Pseudospectral Discretization of Nonlinear Delay Equations: New Prospects for Numerical Bifurcation Analysis. SIAM J. Appl. Dyn. Sys., Vol. 15(1), pp. 1–23, 2016. [3] K.E.M. Church. Validated integration of differential equations with statedependent delay. Commun. Nonlinear Sci. Numer. Simul., 115:106762, 2022. [4] J. Gimeno, J-P. Lessard, J.D. Mireles James and J. Yang, Persistence of Periodic Orbits under State-dependent Delayed Perturbations: Computer-assisted Proofs, SIAM J. Appl. Dyn. Sys., Vol. 22(3), pp. 1743–1779, 2023 [5] A. Gierzkiewicz, R. Szczelina. Sharkovskii theorem for infinite dimensional dynamical systems. Commun. Nonlinear Sci. Numer. Simul., Vol. 146, 108770, 2025. [6] T. Kapela, M. Mrozek, D. Wilczak, and P. Zgliczy´nski. CAPD::DynSys: A flexible C++ toolbox for rigorous numerical analysis of dynamical systems. Commun. Nonlinear Sci. Numer. Simul., 101:105578, 2021. [7] G. Kiss and J.P. Lessard. Computational fixed-point theory for differential delay equations with multiple time lags. J. Diff. Eq., Vol. 252(4), pp. 3093–3115, 2012. [8] J-P. Lessard and J.D. Mireles James. A rigorous implicit C1Chebyshev integrator for delay equations. J. Dyn. Diff. Eq., 33:1959–1988, 2021. [9] M. C. Mackey and L. Glass. Oscillation and chaos in physiological control systems. Science, 197(4300), pp. 287–289, 1977. [10] A. Rauh and E. Auer. Verified Integration of Differential Equations with Discrete Delay. Acta Cybernetica, 25(3), pp. 677–702, 2022. [11] R. Szczelina and P. Zgliczy´nski. High-order Lohner-type algorithm for rigorous computation of Poincar´e maps in systems of delay differential equations with several delays. Found. Comput. Math., 24(4), pp. 1389–1454, 2024.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 90 (a) Bi-directional interconnection. (b) Uni-directional interconnection. (c) Parallel connection. Figure 1: Representation of the structure of interconnected subsystems in CPSs [5]. Contract-Based Controller Tuning In this contribution, we extend the work sketched above by the following multi-stage approach. Firstly, base contracts are derived by means of the procedure published in [5] for an interconnected system with bounded parameter uncertainty and disturbances. Second, we establish an approach for a replacement of controller components in the sense of contract refinement. This approach leads to outcomes of the reachability analysis which are true subsets of the results of the first development stage. As such, the overall structure is still proven to be feasible for the refined system model, however, it also possesses a reduced sensitivity against the considered uncertainties in comparison with the baseline solution. It is further shown how this procedure can be used to validate the use of alternative hardware components (such as actuators or sensors) in a complex interconnected CPS. References [1] J. Alexandre dit Sandretto, A. Chapoutot, and O. Mullier. Formal Verification of Robotic Behaviors in Presence of Bounded Uncertainties. Proc. of 2017 First IEEE International Conference on Robotic Computing (IRC), pages 81–88, 2017. [2] A. Benveniste, B. Caillaud, and R. Passerone. A Generic Model of Contracts for Embedded Systems. 2007. [3] F. Bruns, A. Rauh, and M. Fnadi. Set-Based Assumption-Guarantee Reasoning for Handling Uncertainty in Functionality and Safety Verification of Dynamic Systems. In A. Rauh, B. Finkbeiner, and P. Kr¨oger, editors, Design and Verification of Cyber-Physical Systems: From Theory to Applications, A Springer Nature Computer Science book series, 2025. [4] B. Meyer. Applying ‘Design by Contract’. Computer, 25(10):40–51, 1992. [5] A. Rauh and F. Bruns. A Set-Based Approach to Derive Contracts for Dynamic Behavior in Interconnected Systems, Proc. of the 29th Intl. International Conference on Methods and Models in Automation and Robotics (MMAR), Mi, edzyzdroje, Poland, 2025. Accepted.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 91 Oscillating orbits in the Sitnikov model: equal masses case M. J. Capi´nski1, A. Gierzkiewicz2, and P. Mart´ın3,4 1AGH University of Science and Technology, al. Mickiewicza 30, 30-059 Krak´ow, Poland [email protected] 2Jagiellonian University, ul. Lojasiewicza 6, 30-348 Krak´ow, Poland [email protected] 3Universitat Polit`ecnica de Catalunya, Departament de Matem`atiques & IMTECH, Pau Gargallo 14, Barcelona, Spain 4Centre de Recerca Matem`atica, Campus de Bellaterra, Edifici C, Barcelona, Spain [email protected] Keywords: Sitnikov problem, oscillatory motion, computer-assisted proof. Abstract The Sitnikov three body problem (S3BP) is an example of a spatial elliptic 3BP with oscillatory motion. The configuration has a planar binary consisting of two symmetric bodies of mass mrunning around the common center of mass, and the other body of mass m1, allowed to move along their perpendicular axis of symmetry. Originally [3] m1= 0, which made the system a 1.5 degrees of freedom Hamiltonian problem. Sitnikov showed the existence of oscillatory motion for this case. By ‘oscillatory motion’ we mean an orbit with the mass m1going closer and closer to infinity but always returning to a fixed bounded region. In our case of 0 =m1=m, after rescaling, the system can be seen as a 2 d.o.f. Hamiltonian system with the energy function H(r, ρ, R, y) = 1 2R2+1 2y2+1 ρ2−1 pr2+ρ2−1 4ρ, (1) where ris the position of the third body on the perpendicular axis, ρroughly describes the size of the binary, and R,yare the conjugate momenta. In a joint work with M. Capi´nski and P. Mart´ın, we prove the existence of oscillatory orbits for the Sitnikov 3BP in the case of three equal masses (m1=m). The proof relies on analyzing the stable and unstable invariant manifolds of infinity and their intersections. We construct orbits shadowing these invariant manifolds by the method of correctly aligned windows, which is a modification of the method used in [2]. The proof is computer assisted with the use of CAPD C++ library for rigorous integration of ODEs and computation of Poincar´e maps derived from ODEs [1]. Acknowledgement AG and MC were supported by Polish NCN grant 2021/41/B/ST1/00407.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 92 References [1] T. Kapela, M. Mrozek, D. Wilczak and P. Zgliczy´nski. CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems. Communications in Nonlinear Science and Numerical Simulation, 101:105578, 2021. [2] M. J. Capi´nski, M. Guardia, P. Mart´ın, T. M. Seara, P. Zgliczy´nski et al. Oscillatory motions and parabolic manifolds at infinity in the planar circular restricted three body problem. Journal of Differential Equations, 320:316–370, 2022. [3] K. Sitnikov. The existence of oscillatory motions in the three-body problem. Soviet Physics Doklady, 5:647, January 1961.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 93 Interval Particle Filter for LiDAR-Based Object Tracking Mohamed Fnadi and R´egis Lherbier LISIC – UR4491, Universit´e du Littoral Cˆote d’Opale F-62228, France {mohamed.fnadi,regis.lherbier}@univ-littoral.fr Keywords: Interval Particle Filter, Target Tracking, LiDAR Sensor Introduction Accurate target tracking remains challenging for autonomous systems due to sensor noise and inherent uncertainties. Particle filters (PF) [1, 2] handle non-Gaussian noise effectively but require precise noise and accurate models, while set-membership methods [3, 4] guarantee bounded estimates yet often remain overly conservative. Building on [5], we propose an Interval Particle Filter (IPF) that integrates the probabilistic robustness of PF with the bounded-error guarantees of interval analysis for robust LiDAR-based tracking. The IPF represents the target state using weighted interval particles {[x(i) k], ω(i) k}Np i=1, updated through prediction–correction cycles. As illustrated in Figure 1, the proposed IPF algorithm processes uncertain LiDAR measurements [yk] = {[rk],[αk]}, where [rk] denotes the bounded radial distance to the target’s barycenter and [αk] the orientation deviation between LiDAR and target frames. Predicted box particles at k Predicted box particles at k+ 1 [x(i) k]=[f]([x(i) k−1],[uk−1]),[w(i) k−1]Np i=1 [x(i) k+1]=[f]([x(i) k],[uk]),[w(i) k]Np i=1 X Y xL yL LiDAR sensor x y k x y k+1 Target Trajectory Contraction steps with measurement box [yk]∩[g][x(i) k]Np i=1 [yk+1]∩[g][x(i) k+1]Np i=1 Figure 1: Principle of the IPF filter. During prediction, each particle’s state interval is propagated via the inclusion function [f]([x(i) k−1],[uk−1],[w(i) k−1]), where [uk−1] captures bounded control inputs and [w(i) k−1] represents process noise bounds. In correction, predicted measurements [y(i) k] = [g]([x(i) k]) are compared with actual measurements, and the state intervals are contracted via set inversion as [˜ x(i) k]=[x(i) k]∩[g]−1([yk]⊖[v(i) k]), where [v(i) k] bounds measurement noise. Weights are updated based on the contraction ratio ω(i) k∝ µ([˜ x(i) k])/µ([x(i) k]), and the final state estimate ˆ xkis obtained from the weighted midpoints of contracted particles. Systematic resampling is performed when the effective sample size Neff = 1/P(ω(i) k)2falls below Np/3 to maintain particle diversity.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 94 Main results The IPF was implemented using INTLAB and evaluated against a conventional PF for Autonomous Surface Vehicle (ASV) tracking using LiDAR measurements. Both filters processed identical LiDAR datasets under bounded noise and initial uncertainty. As illustrated in Fig. 2, the IPF achieved superior accuracy (e.g., RMSE of 0.0246 m with 200 particles versus 0.0823 m for the PF with 1000 particles) while providing guaranteed state bounds. However, this came at the cost of higher computation time (e.g., 0.34 s), highlighting the need for further optimization. Figure 2: Comparison of tracking performance: (a) shows the computed trajectories of IPF and PF methods, (b) displays the corresponding lateral errors during ASV tracking. References [1] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp. A tutorial on particle filters for online nonlinear non-Gaussian Bayesian tracking. IEEE Transactions on signal processing, 50(2), pp. 174–188, 2002. [2] B. Fortin, J. C. Noyer, and R. Lherbier. A particle filtering approach for joint vehicular detection and tracking in lidar data. 2012 IEEE Inter Instrumentation and Measurement Technology Conference Proceedings, pp. 391–396, 2012. [3] A. Rauh, M. Lahme, S. Rohou, L. Jaulin, T. N. Dinh, T. Raissi, and M. Fnadi. Offline and Online Use of Interval and Set-Based Approaches for Control and State Estimation: A Review of Methodological Approaches and Their Application. Logical Methods in Computer Science, 2024. [4] M. Fnadi, and A. Rauh. Exponential State Enclosure Techniques for the Implementation of Validated Model Predictive Control. Acta Cybernetica, 26(4), 839–854, 2024. [5] F. Abdallah, A. Gning, and P. Bonnifait. Box particle filtering for nonlinear state estimation using interval analysis. Automatica, 44(3), pp. 807–815, 2008.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 95 Friday, September 26, 2025 A01-0-006 11:00–12:30 Regular Session B: Optimization 11:00–11:30 Mih´aly Gencsi and Bogl´arka G.-T´oth: Improvements of the Geometrical Test in Interval Branch and Bound methods 11:30–12:00 Verlein Radwan, Simon Rohou and Gilles Trombettoni: Exhaustive Interval-based 2-D Shape Registration Under Similarity Transformation 12:00–12:30 Ma¨el Godard, Luc Jaulin and Damien Mass´e: Adaptative parallelepipedic approximation of the image of a set by a nonlinear function 14:00–15:30 Regular Session B: Optimization 14:00–14:30 Yuki Uchino, Katsuhisa Ozaki and Toshiyuki Imamura: High-Performance Emulation of Matrix Multiplication using INT8 Matrix Engines and its Error Analysis 14:30–15:00 Lorenz Gillner and Ekaterina Auer: Efficient Acceleration Strategies for Interval Branch-and-Bound Type Methods 15:00–15:30 Diego Romano, Ekaterina Auer, Francesco Gregoretti and Lorenz Gillner: GPU-Accelerated Algorithmic Differentiation For Reliable Computing: Comparing Different Architectures
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 96 Improvements of the Geometrical Test in Interval Branch and Bound methods Mih´aly Gencsi, Bogl´arka G.-T´oth University of Szeged Department of Computational Optimization H-6720 Szeged, Hungary {gencsi,boglarka}@inf.szte.hu Keywords: Interval Methods, Branch and Bound method, Fritz-John optimality conditions, Geometrical Test Introduction In many real-world applications, finding a guaranteed solution is crucial. We focus on the following nonlinear inequality-constrained global optimization problem, minimize x∈y⊆Rnf(x) subject to gi(x)≤0, i ∈Mc, (1) where f, gi:Rn→R, i ∈Mcare continuously differentiable nonlinear functions. The interval box y= [y, y] defines the general bound constraints. These constraints can be expressed as piu(x) = xi−yiand pil(x) = yi−xi, which can be written compactly as pj(x)≤0, j ∈Mb, where Mb={iu, il|i= 1, . . . , n}. In this study, we improve the efficiency of the IBB method using optimality conditions. We improve the Advanced Geometrical Test to discard non-optimal boxes and avoid using optimality conditions when they are ineffective. The Normalized Interval Fritz-John Condition System The normalized interval Fritz-John condition system (NIFJ-CS), as used in [1], defines a system of interval linear equations for a given box x: ϕ(t) = µ0+X i∈B∪C µi−1 µ0·∇f(x) + Pi∈Bµi·∇pi(x) + Pj∈Cµj·∇gj(x) µi·pi(x)i∈B µj·gj(x)j∈C = 0,(2) where f(x),pi(x), gj(x) are the inclusion functions and ∇f(x), ∇pi(x), ∇gj(x) are the inclusions of the gradients. The system includes only active constraints (B= {i∈MB|0∈pi(x)},C={i∈MC|0∈gj(x)}). The variables are t= [x,µ]T. Several approaches to solving Fritz-John optimality conditions have been proposed in the literature. As discussed in [2, 3], two main methods are: computing only the bounds for the Lagrange multipliers (Lagrange estimator method) and enclosing both the multipliers and the optimal point location (Lagrange estimator+NIFJ-CS method). Both methods are computationally expensive and often ineffective. Therefore, an efficient modification is needed.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 97 Advanced Geometrical Test Instead of applying the Lagrange estimator or combining it with the NIFJ-CS method, we perform a preliminary test, the Advanced Geometrical Test. This test discards the nonoptimal boxes or avoids the unnecessary use of the optimality conditions. The Advanced Geometrical Test is based on the geometrical interpretation of the optimality conditions [4]. Geometrically, in the interval settings, the optimality conditions require that there exists a direction in the negative cone of the objective’s gradient, CF ={d∈Rn|d∈ −µ0∇f(x), µ0≥0}, which lies in the conic hull of the gradient enclosures, CH =(d∈Rn|d∈X i∈B µi∇pi(x) + X i∈C µi∇gi(x), µ ≥0, µ = 0). So, if CF ∩CH =∅, then the interval box xcan contain a local optimizer. The test performs step-by-step checks on the conic hull CH and the cone CF, which allows for early termination and saves computation time. The procedure begins with a check to see if the gradient, denoted by ∇f(x), contains zero, or if any active constraint’s gradient has zero in its interior. If either condition is satisfied, the box cannot be reduced. Next, we analyze the sign patterns of the gradient enclosures in each dimension and classify them as “+”, “–”, “0+”, “0–”, or “+–”. If the signs of the conic hull and the cone differ in any dimension, the necessary optimality conditions cannot hold, and the box is discarded. Then, we check if the gradient enclosures of the active constraints cover all orthants. In this case, CH is full, so the box cannot be reduced. For an active constraint, we compute the inclusion of the Lagrange multipliers in each dimension. If the intersection of these inclusions is empty, the box is discarded. For multiple active constraints, we compute the same inclusions using the interval hull of the active gradients. Again, if the intersection is empty, the box is discarded. Lastly, we consider pairs of dimensions where at least one has sign-consistent gradient enclosure. For each pair, we compute the slope intervals Mf and Mg. If Mf∩Mg=∅, the box is eliminated. Efficient pairing strategies reduce the running time of the pairing method. We show that the above Advanced Geometrical Test is very efficient on a large benchmark. The best variant can save more than 40% of the computational time on average using the designed test. References [1] M. Gencsi and B. G.-T´oth. The Fritz-John condition system in interval branch and bound method. Annales Mathematicae et Informaticae, 58:56–68, 2023. [2] E. Hansen and G.W. Walster. Global Optimization Using Interval Analysis: Revised And Expanded. CRC Press, 2003. [3] E.R. Hansen and G.W. Walster. Bounds for Lagrange multipliers and optimal points. Computers & Mathematics with Applications, 25(10):59–69, 1993. [4] M. Gencsi and B. G.-T´oth. Efficient use of optimality conditions in Interval Branch and Bound methods. EURO Journal on Computational Optimization, 13, 2025.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 98 Exhaustive Interval-based 2-D Shape Registration Under Similarity Transformation Verlein Radwan1, Simon Rohou2and Gilles Trombettoni1 1LIRMM - Lab. d’Informatique, de Robotique et de Micro´electronique de Montpellier 161 Rue Ada, 34095 Montpellier, France {vradwan,gilles.trombettoni}@lirmm.fr 2ENSTA, Institut Polytechnique de Paris, Lab-STICC 2 Rue Fran¸cois Verny, 29200 Brest, France [email protected] Keywords: Procrustes Analysis, Registration, Separators, Interval Analysis Motivation The registration of two sets aims at getting all the transformation parameters that map them. This work covers the case of bounded subsets of R2and similarity transforms, consisting of a composition of translation, rotation and uniform scaling. The classical Procrustes-like approach [1] to this problem is to switch from the original 4D transformation to a 1D orientation problem. Since rotational symmetries induce several solutions for the orientation problem, state-of-the-art methods based on local optimization are de facto unsuitable when they appear. We propose a set-membership approach capable of approximating all solutions by reproducing the Procrustes-like approach and solving the orientation problem, thanks to set manipulation using descriptive operators called separators [2]. A Procrustes set-membership approach The algorithmic sequence used to perform the registration of two sets Aand Bstarts by centering and normalization steps, followed by rotational mapping. The usual centers and scaling factor come from the fact that Aand Bare represented as point clouds in state-of-the-art methods. In our set-membership approach, we have to find appropriate substitute for these (centers and scaling factors) parameters. Proposition 1. The minimum enclosing circle [3] provides different but suitable parameters for the centering and normalization. Thus, the whole algorithm consists of the following steps: 1. Find the minimum circles centers cA 1,2,cB 1,2and radii cA 3, cB 3. (Fig. 1.A), 2. Describe the normalized sets AN:= A−cA 1,2/cA 3,BN:= B−cB 1,2/cB 3. (Fig. 1.B–C), 3. Describe the set of possible rotations Θ = {θ|BN= Rθ·AN}, where Rθis the rotation matrix of parameter θ. (Fig. 1.D).
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 99 A B C D Figure 1: Algorithm in 4 steps: A: Identification of smallest circles (centers and sizes). B: Centering of sets. C: Normalization of sets. D: Identification of rotation parameters. Approximation using separators Usually involved in a branch and separate algorithm, a separator [2] is an algorithmic operator capable of representing a set: from an initial domain, it can remove nonsolution parts (contraction) as well as parts containing only solutions. Using the Codac library (www.codac.io, [4]) which offers a catalog of separators, it is possible to define a new separator as a sequence of separators handling constraints involving sets and set operations. To implement the Procrustes approach, we first give the system of constraints that defines the problem, and for each variables we propose a separator consistent with it. As a result, the algorithm returns separators describing respectively all possible sets Θ,K,Tof rotations θ, (unique) uniform scaling kand translations t. Finally, by slightly modifying the chain of separators, we will see that the same approach allows similarity transformations to be approximated completely, including reflection symmetries. References [1] P.H. Sch¨onemann and R.M. Carroll. Fitting one matrix to another under choice of a central dilation and a rigid motion. Psychometrika, vol 35, 245–255, 1970. [2] L. Jaulin and B. Desrochers. Introduction to the Algebra of Separators with Application to Path Planning. Engineering Applications of Artificial Intelligence, vol 33, 141–147, 2014. [3] J.J. Sylvester. A question in the geometry of situation. Quarterly Journal of Pure and Applied Mathematics, vol 1, 79–80, 1857. [4] S. Rohou and B. Desrochers and F. Le Bars. The Codac Library. Acta Cybernetica, Special Issue of SWIM 2022, vol 26, 871–887, 2024.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 106 of GPUs into existing CPU-based code. Additionally, we assess our findings using a close-to-life application. Finally, we consider the implications of GPU acceleration with respect to power consumption. Although GPUs offer substantial performance improvements with each new generation, they also require significant amounts of energy. To quantify our energy-related costs, we not only benchmark throughput performance but also perform detailed energy measurements for both CPU-only and GPU-accelerated interval implementations. Unlike our previous studies [1, 6], which relied on software-based estimates, this work presents our results using hardware-based measurements, offering a more accurate assessment of the trade-offs between performance and power consumption in parallelized IBB methods. References [1] E. Auer, A. Ahrens and L. Gillner. GPU-Based Interval Optimization in the Context of Optical MIMO Systems. In: R. Wyrzykowski, J. Dongarra, E. Deelman, and K. Karczewski (Eds.), Parallel Processing and Applied Mathematics, 360–374, 2025. [2] E. Auer, A. Rauh and L. Gillner. Parameter Identification for Cooperative SOFC Models on the GPU. In: T. N. Dinh, A. Rauh, S. Z. Yong, Z. Wang (Eds.), Set-Valued Approaches to Control and Estimation of Uncertain Systems – Theory and Applications, to be published. [3] Z. Bag´oczki and B. B´anhelyi. A Parallel Interval Arithmetic-based Reliable Computing Method on a GPU. Acta Cybernetica, 23:491–501, 2017. [4] P.-D. Beck and M. Nehmeier. Parallel Interval Newton Method on CUDA. Applied Parallel and Scientific Computing, 454–464, 2013. [5] C. Collange, M. Daumas, D. Defour. Interval Arithmetic in CUDA. In: W. W. Hwu (Ed.): GPU Computing Gems Jade Edition, 99–107, 2012. [6] L. Gillner and E. Auer. GPU-Accelerated, Interval-Based Parameter Identification Methods Illustrated Using the Two-Compartment Problem. Acta Cybernetica, 26:913–932, 2024. [7] E. Hansen and G. W. Walster (Eds.). Global Optimization Using Interval Analysis: Revised and Expanded. CRC Press, 2004. [8] L. Jaulin and E. Walter. Set Inversion via Interval Analysis for Nonlinear Bounded-error Estimation. Automatica, 29:1053–1064, 1993. [9] R. E. Moore, R. B. Kearfott, and M. J. Cloud. Introduction to Interval Analysis. Society for Industrial and Applied Mathematics, 2009. [10] D. P. Sanders and V. Churavy. Branch-and-bound interval methods and constraint propagation on the GPU using Julia. SCAN, 2020.
Efficient Acceleration Strategies for Interval Branch-and-Bound Type Methods Lorenz Gillner and Ekaterina Auer University of Applied Sciences Wismar Department of Electrical Engineering and Computer Science Motivation Classical interval algorithms for problems such as root-finding, global optimization and parameter estimation iteratively apply branching and bounding/pruning on an initial search domain until a satisfactory result is reached. These types of algorithms are often restricted to problems of low dimensions due to their dependence on space subdivision. As a result of the growing interest in GPU applications for scientific computing in general, the topic of GPU-accelerated interval methods has also gained more attention in recent years. We explore ways to enhance the speed of interval branch-and-bound type methods using GPUs to scale them beyond their current limits. Given today’s growing demand for computing power, we must evaluate accelerated algorithms based not only on their speed, but also on their energy requirements. GPU Computing Key considerations: nFit tasks to thread block model nReduce main memory access nMaximize occupancy åAvoid idle threads nKernels operate in warps (groups of 32 threads) Thread Block Model Grid Block Thread Warp Parallelization of Branching Algorithms Two main categories of parallelization for tree-like procedures: Node Level Work distribution across a coarse grid (subtasks) Leaf Level Subdivision of a problem into fine grid (brute-force) General search space reduction scheme used previously [1, 2, 3]: 12 GPU 34 Drawbacks: High memory saturation; scalability limited by problem dimensions Improvement: Unify steps 1 and 2 to maximize computational efficiency if f([x]) ⊂[y]:keep [x] elif f([x]) ∩[y] = ∅:reject [x] else:bisect [x] SIVIA (Set Inversion Via Interval Analysis [4]): Good example for branching-based interval algorithms suitable for parallelization +We consider 3 basic parallelized versions, as well as combinations thereof: NSIVIA Multiple (independent) instances of SIVIA in parallel over a coarse grid PSIVIA Massively parallel, one-time evaluation of single boxes over a fine grid VSIVIA Vectorized »bisect – evaluate – partition« cycle [5] Boost.Interval: New GPU Mode To streamline development, we added GPU support to the Boost.Interval library: äDirect rounding mode for GPU arithmetic ßfast operations äExtensive use of generic programming ßsame source code for CPU and GPU äFirst cross-platform, GPU-compatible interval library for CUDA C++ /* Classic interval type definition */ typedef interval<double, policies<save_state<rounded_transc_std<double> >, checking_base<double>>>host_interval; /* Interval data type definition for the GPU */ typedef interval<double, policies<save_state_nothing<rounded_transc_gpu<double> >, checking_base_gpu<double>>>device_interval; /* Cross-platform interval type (without bounds checking) */ typedef interval<double, policies<save_state_nothing<rounded_transc_exact<double> >, checking_base<double>>>cross_interval; Githubhttps://github.com/lorenzgillner/boost-interval Performance Benchmarks Benchmarks: SIVIA for the Griewank test function g(x) = 1 4000 n X i=1 x2 i− n Y i=1 cos xi √i!+ 1,n∈ {1, . . . , 4} nOptimization of the method itself åAll results are qualitative the same n16 variations of SIVIA implementations nComparison of parallelization frameworks åOpenMP, OpenMPI, TBB and CUDA nTests on consumer grade hardware Power Profiling äPer-component measurements of power consumption during benchmarks 20 40 60 80 Time (s) 0 50 100 150 200 Power (W) CPU GPU Sync Results 12 3 4 1 10 100 Dimension Speedup NSIVIA PSIVIA VSIVIA äChoice of parallelization method depends on the problem dimensions äBrute-force methods are best suited for small problems äA-priori subdivision combined with vectorization yields best time to solution ä675speedup on a four-dimensional problem compared to traditional version åReduction of energy requirements by 90% References: [1] L. Gillner and E. Auer. GPU-accelerated, interval-based parameter identification methods illustrated using the two-compartment problem. Acta Cybernetica, 2024. [2] M. B. Eriksen and S. Rasmussen. GPU accelerated parameter estimation by global optimization using interval analysis, 2013. [3] K. Nasiotis. Accelerating SIVIA (set inversion via interval analysis). an interval set membership technique to evaluate the generalization of neural classifiers, 2023. [4] L. Jaulin and E. Walter. Guaranteed nonlinear parameter estimation from bounded-error data via interval analysis. Mathematics and computers in simulation, 1993. [5] S. Tzeng and J. D. Owens. Finding convex hulls using quickhull on the GPU. 2012. 20th International Symposium on Scientific Computing, Computer Arithmetic, and Verified Numerical Computations, Sept. 22–26, 2025, Oldenburg, Germany
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 108 GPU-Accelerated Algorithmic Differentiation For Reliable Computing: Comparing Different Architectures Diego Romano1, Ekaterina Auer2, Francesco Gregoretti1and Lorenz Gillner2 1Institute for High Performance Computing and Networking National Research Council, Italy {diego.romano, francesco.gregoretti}@cnr.it 2Department of Electrical Engineering University of Applied Sciences Wismar, D-23966 Wismar, Germany {ekaterina.auer,lorenz.gillner}@hs-wismar.de Keywords: Interval Methods, Algorithmic Differentiation, GPU As the need for trustworthy simulations grows, reliable computations are becoming increasingly important in such varied areas as robotics, medical imagining, computational biology, and many others. In such applications, accurately propagating bounded epistemic uncertainty is crucial for achieving dependable results in both static and dynamic systems. Beyond simulation, parameter identification under uncertain conditions plays a vital role in solving real-world engineering and scientific problems. The demand for reliability is offset by an equally important requirement of efficiency, often resulting in a compromise that prioritizes one over the other. Interval methods [3] have emerged as a valuable tool in the context of reliability, since they inherently propagate bounded uncertainty while fulfilling their primary purpose of result verification. This kind of verification involves providing a guaranteed enclosure of the exact result, despite potential rounding, discretization, or method errors that may affect computer-based calculations. An important step in the verification of higher-level numerical methods, such as those used for solving differential equations or global optimization, is the computation of exact partial derivatives of the underlying functions. While users can provide the code for derivatives manually or with the aid of a computer algebra system, it is more desirable to compute derivatives automatically within the normal program code. Algorithmic differentiation [2] offers a well-established approach to achieve this goal. Over the past three decades, numerous tools have been developed to support its implementation, primarily targeting serial execution on the CPU. While algorithmic differentiation simplifies the process from the point of view of human effort, it must not necessarily result in improved overall runtime performance of a simulation. Parallelization has been successfully employed within the high-performance computing framework to accelerate traditional result verification methods. The use of graphic processing units (GPUs) can be particularly beneficial, owing to their high computational throughput and relatively low cost. Nevertheless, the GPU architecture diverges significantly from its CPU counterpart, necessitating innovative approaches to porting CPU code to the GPU in order to fully exploit the GPUs’ capabilities.
SCAN 2025, Sept. 22-26, 2025, Oldenburg, Germany 109 This poses a particular challenge for algorithmic differentiation, as only a few GPU implementations currently exist1, and none are designed to support interval data types. To address this gap, we have recently ported the popular template-based CPU tool FADBAD++ to the GPU and implemented a differentiation-algebra-based approach for the first and second derivative2. Although FADBAD++ has the big advantage of enabling the computation of derivatives of any order on the CPU, a naive porting to the GPU proves to be inefficient [1]. In this contribution, we initiate a systematic approach to optimizing the implementation of algorithmic differentiation on the GPU. By analyzing the performance characteristics of various GPU architectures, we aim to pinpoint the primary bottlenecks and identify opportunities for improvement. To this end, we examine the computational efficiency of a classic global optimization problem using the Griewank function, alongside a real-world case from the field of multiple-input multiple-output (MIMO) systems, a key technology in broadband wireless communications. References [1] E. Auer, A. Ahrens, and L. Gillner. GPU-based interval optimization in the context of optical MIMO systems. In: R. Wyrzykowski, J. Dongarra, E. Deelman, and K. Karczewski (Eds.), Proceedings of Parallel Processing and Applied Mathematics, 2025. [2] A. Griewank. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. Frontiers in Applied Mathematics Series, 2000. [3] R.E. Moore, R.B. Kearfott, and M.J. Cloud. Introduction to Interval Analysis. Society for Industrial and Applied Mathematics, 2009. 1e.g., https://github.com/vgvassilev/clad 2https://github.com/lorenzgillner/tada