Advantages of the usage of the Infinity Computer for reducing the Zeno behavior in hybrid system models

To capture the dynamics of modern Cyber-Physical Systems, hybrid system models are introduced to combine their continuous dynamics with the discrete ones. Unfortunately, one important negative issue can affect hybrid system models: the so-called Zeno phenomenon, which results in an infinite number of discrete transitions in a finite amount of time occurring during the model’s simulation that leads to inconsistent results. In this context, the paper investigates the use of a recently proposed numerical algorithm, based on the Infinity Computer methodology, to handle the Zeno phenomenon and evaluate it with respect to standard numerical methods by considering the hybrid system models of two exemplary Cyber-Physical Systems: the Water tanks and the Thermostat.


Introduction
To design and analyze Cyber-Physical Systems (CPSs), presenting both continuous and discrete dynamics, hybrid system models have been introduced and are widely adopted [see, e.g., Bocciarelli et al. (2019), , Grossman et al. (1993); Lunze and Lamnabhi-Lagarrigue (2009)]. A hybrid system can be defined as a system whose dynamic is regulated through continuous and discrete behaviors. Ordinary differential equations (ODEs) are often used to govern the continuous behavior, whereas control graphs can be used to describe the discrete behavior. The values of continuous variables in a particular discrete mode define the state of a hybrid system [see, e.g., Bouskela et al. (2021), Platzer (2008)]. More in detail, continuous states can change their values in two ways: (i) through a discrete transition or (ii) based on the differential equations results. Discrete states can only alter values through discrete transitions. As a consequence, hybrid systems can be considered an efficient formalism for modeling and simulating CPSs, where digital components interact with a physical environment, and their dynamics are regulated through continuous and discrete behaviors.
One of the most important issues unique to hybrid systems is the Zeno phenomenon that can happen when the system undergoes an unbounded number of discrete transitions in a finite and bounded length of time [see, e.g., Johansson et al. (1999), Zhang et al. (2001)]. To explain this phenomenon and its effects, suppose to model and simulate a ball that is dropped from a predefined height h 0 > 0 with an initial velocity v 0 = 0. The ball in free fall hits the ground after a certain time t in partially elastic way, loses energy, bounces back into the air, and then starts falling again until it stops. Starting from this physical system, a hybrid automata is designed in order to model the behavior of the bouncing ball [see Falcone et al. (2022), Johansson et al. (1999) for details]; in particular, ODEs are used to model the variation of h and v during the motion of the ball, whereas transitions rules are introduced to model the hit of the ball with the ground causing the ball to bounce. As the ball loses energy when it hits the ground, with the progress of time, a large number of discrete transitions occur progressively in the hybrid automata in smaller and smaller time intervals; thus, the automata experiences the Zeno phenomenon. This behavior leads to have, in a simulated environment, incorrect data (e.g., negative values for h) resulting in unfeasible behavior (e.g., the ball goes below the ground level) [see Fritzson (2014) for details].
It is important to note that the real-world CPSs do not show Zeno behavior, but their models may exhibit Zeno executions due to modeling abstractions [see, e.g., Heymann et al. (2005)]. The Zeno phenomenon is difficult to characterize and mitigate because it leads to incorrect simulation results since the system behavior is intrinsically ill-conditioned beyond the point in time (called Zeno point) in which such a phenomenon starts [see, e.g., Branicky (2005), Zhang et al. (2000)]. In this situation, it is necessary to carefully check the prerequisites for discrete events' transitions.
Different approaches and techniques for addressing the Zeno phenomenon have been proposed in the literature. In order to assess the effect on real-world systems, some researchers have examined this issue when hybrid system automata are derived as abstractions of the underlying physical systems and defined the rules associated with this phenomenon [see, e.g., Ames and Sastry (2004), Falcone and Garro (2020), Heymann et al. (2005), Lygeros et al. (2003)]. Other researchers have used regularization and sliding mode approaches to extend the simulation of Zeno hybrid systems beyond the Zeno point. For example, in Johansson et al. (1999), the authors delineated the characteristics of hybrid automata and employed regularization approaches to simulate and evaluate these automata beyond the Zeno point. However, in order to apply regularization approaches it is necessary to alter the original system by perturbing it with a small quantity to obtain a non-Zen solution. If = 0, the regularized system tends to be the original one. However, such perturbation may invalidate the notion of instantaneous discrete transitions and the simulation performance may degrade. The other approach to handle the Zeno phenomenon regards the use of a sliding mode algorithm that is based on detecting regions on the switching manifold on which the Zeno phenomenon occurs, and then forcing the system to slide on the manifold in these regions [see, e.g., Utkin (2013), Weiss et al. (2015), Yu et al. (2011)]. When infinitely fast transitions occur, a smooth sliding action takes place on the switching surface to eliminate the Zeno phenomenon. It is also possible to modify the original hybrid system by adding the extra mode to represent the dynamics during the mode transitions, but sliding mode algorithms can only be applied to study special classes of hybrid systems [see, e.g., Biák et al. (2013), Filippov (2013]. From a numerical point of view, extending the simulation of Zeno hybrid systems beyond the Zeno point means the definition of fast and efficient procedures for the determination of zero-crossings of the transition conditions. Indeed, classical numerical algorithms for simulating hybrid system models generate observations t k = t k−1 + t k−1 , k ∈ {1, 2, ...}, of time t ∈ [t 0 , T ] in different ways. The standard method consists of using a fixed stepsize t k = t = ω with k ∈ {0, 1, 2, ...}, but there exist other methods for generating observations dynamically, such as MATLAB/Simulink procedures using Dormand-Prince methods [see, e.g., Kimura (2009)]. However, the main objective of these methods is to solve the ODEs used to describe the system in an efficient way, and not to determine the zero-crossings. This means that these methods are useful to simulate hybrid system models under a fixed state and a fixed set of ODEs, but they do not allow determining whether a zero-crossing occurs between two consecutive observations t k and t k+1 .
To overcome the issues presented above, in Falcone et al. (2022), an efficient method, called Infinity Computer algorithm using Runge-Kutta method of the fourth order (IC_RK4), to simulate hybrid system models by generating observations t k , k ∈ {0, 1, 2, ...}, dynamically has been proposed. The proposed method allows studying better the regions, where zero-crossings can occur, without spending a lot of computational resources on "stable" regions. This is done because, in a hybrid system, the continuous behavior is performed as long as invariants remain, whereas discrete transitions (or events) occur when a specific jump condition is satisfied. As a consequence, in the IC_RK4 method, a jump condition is associated with zero-crossings of a given function g(x(t)). This indicates that an event occurs at a time t, when there have been generated the system variables x = x(t), such that g(x(t)) = 0 [see, e.g., Ames et al. (2006)]. The IC_RK4 method has been developed taking into consideration the following characteristics: -Zero-crossing detection. Zero-crossings are detected by evaluating approximations of the original function g(x(t)), without resorting to the original hybrid system. This characteristic allows efficient detection of zero-crossings without requiring additional computing resources, even if the original system is difficult to simulate; -Dynamic generation of observations. The built-in algorithm allows generating observations more frequently around zero-crossing regions, and less frequently elsewhere, avoiding unnecessary computations in "stable" regions, where no zero-crossings can happen; -Extensibility. It can be easily extended to integrate any fast and efficient method for detecting zero-crossings, such as global optimization algorithms [see, e.g., Casado et al. (2002), Molinaro and Sergeyev (2001)].
The IC_RK4 method adopts the Infinity Computer that represents a new kind of supercomputer that allows one to work numerically with infinite and infinitesimal numbers [see, e.g., Sergeyev (2010Sergeyev ( , 2017] in a novel framework different w.r.t. nonstandard analysis [see Sergeyev (2019)].
The Infinity Computer allows one to work numerically with finite, infinite, and infinitesimal numbers. It introduces a new numeral ①, called grossone, that represents the number of elements of the set, N, of natural numbers and it is used as the radix of a positional numeral system, in which a number C (finite, infinite, or infinitesimal), known as grossnumber, is expressed in the following form [see, e.g., Falcone et al. (2020b), Sergeyev (2017 for details]: where c i , i = 1, ..., N , named grossdigits, are finite positive or negative real numbers, whereas p i , i = 1, ..., N , known as grosspowers, are arranged in decreasing order; thus, p 1 > p 2 > · · · > p N , and they can in turn be finite, infinite, or infinitesimal quantities of the same form (1). According to this methodology, one of the following three situations can happen: p 1 > 0. The number C is infinite; thus, ① = 1① 1 is the basic infinite number; p 1 < 0. The number C is infinitesimal; thus, ① −1 represents the simplest infinitesimal; p 1 = 0. The number C is finite. Specifically, if N = 1, then C is known as purely finite, i.e., there are no infinite nor infinitesimal quantities.
Arithmetic operations between grossnumbers are carried out in the Infinity Computer as follows. Given the grossnumbers A, B, and C specified as: The addition operation is carried out by including in the result The operation of subtraction is directly derived from addition.
The result of multiplying A and B produces a grossnumber C, defined as follows: where Finally, the division operation of A by B yields a result C with a reminder R, where the first grossdigit is c k K = a l L /b m M , the highest exponent is k K = l L − m M , and the first partial reminder R * is produced as follows: If R * = 0 or the default accuracy is achieved, the division operation is complete; otherwise, the computation is restarted by substituting A with R * .
In this paper, the IC_RK4 method proposed in Falcone et al. (2022) is further exploited to study complex Zeno hybrid systems. More specifically, two well-known Zeno hybrid systems, i.e., Water tanks [see Johansson et al. (1999)] and Thermostat [see Johnson et al. (2004)] are considered in this work. Both the systems are modeled with constant and nonlinear functions in order to stress the IC_RK4 method and evaluate its performance in detecting zero-crossings through a dynamical generation of time observations by using the infinite quantity ① [see Sergeyev (2017) for its detailed description]. To show the validity of the IC_RK4 method, these hybrid systems have been studied and results are gathered from simulations and compared with the standard method.
The rest of the paper is structured as follows. Section 2 provides an introduction to hybrid systems and the IC_RK4 method. Section 3 presents the Water tanks and Thermostat hybrid systems taken from the literature that exhibit Zeno behavior. Section 4 describes the conducted numerical experiments on such systems along with a comparison of simulation results obtained with the exploited IC_RK4 method and a standard one. Conclusions are presented in Sect. 5. Finally, some additional results related to the conducted experiments are reported in Appendix.

Hybrid system models and the IC_RK4 method
The paper uses notions and concepts from hybrid systems and the IC_RK4 method along with the related algorithm and concepts, as presented in the following subsections.

Hybrid system models
The structure of a hybrid system model can be formally represented as [see Lunze and Lamnabhi-Lagarrigue (2009)] follows where X = R n is the continuous state space; Q is a finite set of discrete states; f is a set of vector fields describing the continuous dynamics for all q ∈ Q; Z is a set of initial states; δ is the discrete state transition function; K is a set of guards describing when a discrete state transition occurs; R is a reset map defining the state jumps; I invariants of the discrete states.
The discrete behavior is modeled through a control graph G = (Q, E), where the set of vertices Q represents the discrete states (also called "operation modes"), whereas the .., n]} represent state transitions managed by the function δ : Q × R n → Q, which determines the discrete successor state q j when the system is in the discrete state q i . Every discrete state q i ∈ Q is linked to a vector field f : Q × R n → R n that belongs to it. The evolution of the continuous state in the discrete state q(t) is described by the equationẋ(t) = f (q(t), x(t)).
The collection of initial states Z ⊂ Q × R n is defined within the hybrid system model. Through the associated vector field f (q, ·), a discrete state q affects the continuous dynamics, while the collection of guards K serves as a representation of the impact of continuous dynamics on the evolution of discrete states. Thus, given a continuous state space X and K ⊂ X , a guard k ∈ K is a region where a discrete state transition may take place if the state x is in k.
The reset map R : Q × Q → 2 R n × 2 R n and the invariants of the discrete states I : Q → 2 R n allow to include state jumps and complete the representation of the interaction between the continuous and the discrete dynamics. Specifically, each discrete state q has associated an invariant i ∈ I that describes the conditions that the continuous state x has to satisfy at q. Invariants and guards play complementary roles; specifically, invariants characterize when a transition must take place, whereas guards represent "enabling conditions" that identify when a particular transition may take place. R is a function that specifies how new continuous states are related to previous ones for a given transition.

The IC_RK4 method and related modeling and simulation process
To design and simulate dynamic systems that exhibit the Zeno phenomenon, a modeling and simulation (M&S) process has been defined. This M&S process consists of a set of interconnected activities for the management of such systems. Starting from the description and requirements of a Zeno dynamic system (input), this M&S process allows one to produce an enhanced version of the system that can be simulated by managing the Zeno phenomenon (output). Figure  The M&S process, of which the IC_RK4 method is part, involves two experienced engineers: system engineers and modeling and simulation engineers. System engineers are in charge of delineating the requirements and characteristics of the real system, in a semi-structured way, by using specialized software and technical documents to define the structure and behavior of the system along with its components, whereas modeling and simulation engineers, starting from the technical reports and resources produced by system engineers, define the corresponding hybrid system automata, simulate it, and then gather the simulation results. The results coming from the performed experiments are jointly evaluated by both engineer experts to evaluate the system performance and explore the system's design before building it.
The M&S process begins inside the swimlane of the system engineers, where the real-world dynamic system is defined by capturing the system's requirements and definitions, and the stakeholders' needs in technical documents ("Dynamic System Requirements" and "Dynamic System Description"). This step ensures that the dynamic system's structure and behavior are proper, rational, and effective by making it explicit and complete in its components.
The obtained technical documents represent the input of the "Requirements Representation" and "Dynamic System Components" steps. In the first step, the system's requirements are formally defined to automate their verification through simulation, whereas the components that make up the system are specified in the "Dynamic System Components" step. After defining the hybrid system in terms of requirements and components, the whole system architecture and behavior are derived in the "Dynamic System Architecture and Behavior" step, and the system specifications are stored into a repository. At these steps, two languages are adopted to formalize the dynamic system: the Unified Modeling Language (UML) and/or Systems Modeling Language (SysML) (Johnson et al. 2007).
UML is a modeling language developed and maintained by the Object Management Group (OMG) that allows one to specify, design, implement, and document the artifacts of software systems. UML is not a programming language but it supports the entire life cycle of a software system in different application domains (e.g., finance, aerospace, and industry 4.0). UML provides visual notations, based on diagrams, that are not only intended for developers but also for stakeholders, and anybody involved in the project. The provided diagrams allow one to capture the characteristics of a system along with the relations among its components. SysML is a general-purpose architecture modeling language, maintained by OMG, for system engineering applications. The SysML language supports the definition, analysis, design, verification and validation of a wide range of systems and systems of systems. The SysML language allows system engineers to create effective SysML models able to capture the characteristics of a system at different levels of abstraction. SysML extends a subset of the UML diagrams by using the UML profile mechanism and provides additional diagrams for managing requirements and parametric constraints. As a result, SysML is more adaptable and expressive than UML for modeling dynamic systems.
Upon obtaining the formal representation of the hybrid system automata H S in terms of UML/SysML diagrams, control graph, and differential equations, the "g(x) function definition" step is executed by modeling and simulation engineers inside the corresponding swimlane. In this step, a g(x) function that captures the discrete transitions of H S is defined and included in the hybrid system automata. This step produces as output an extended version of H S, i.e., H S that is evaluated by system engineers for approval. If the design is approved, the "Dynamic System Simulation" step is performed by taking two input arguments: (i) a property file, Algorithm 1: Simulation of a hybrid system automata i.e., "Properties" that delineates the parameters required to configure the simulation scenario (e.g., simulation time step, solver type, log folder, etc.) and (ii) the extended hybrid system automata H S . Inside this step, the simulation is performed according to Algorithm 1. Finally, in the "Simulation Result Evaluation" step, the simulation results are analyzed to study the system's behavior and evaluate different design alternatives before building it.
Algorithm 1 allows one to execute the simulation of a hybrid system automata H S also in the presence of the Zeno behavior. The algorithm takes as parameters: (i) the infinitesimal quantity ① −1 provided by the Infinity Computer [see Sergeyev (2010) where t min represents the smallest time value between two discrete transitions, and t max is the time interval between two consecutive observations on "good" and "stable" regions; (iii) the simulation stop time T ; and, finally, (iv) the initial discrete state of H S , q.
Until the simulation of H S is completed, the approximation is calculated by using the infinitesimal quantity ① −1 (see Algorithm 1-line 5). Note that x is an n-dimensional vector, whereas ① −1 is an infinitesimal scalar. Thus, all the arithmetic operations are componentwise vector operations, and they are performed numerically by adopting the Infinity Computer and not in a symbolic way [as it is done in standard methodologies, e.g., in Shamseddine and Berz (2000)].
Once the value x k+1,① −1 is derived, the approximation g k+1,① −1 is calculated, and then, the function g(y), y ∈ [0, t], is defined using the coefficients of ① −1 extracted from g k+1,① −1 (see Algorithm 1-lines 6, 7): After the search for zero-crossings of g(y) for y ∈ [0, t], two situations can arise: no zero-crossings No zero-crossings have been detected, this means that the function g(y) is positive over the whole interval [0, t]. In this situation, the interval [t k , t k + t] is considered "stable"; thus, it is no necessity to use a smaller stepsize t to better analyze this interval. The value x k+1 is obtained by using the function FinitePart(t), which is provided by SSIC, at the point t = x k+1,① −1 . Since the interval is "stable" for the next simulation step, the value of t is set to t max (see Algorithm 1-lines 9-12); zero-crossings detected Zero-crossings have been detected, this means that the function g(y) can be negative in the interval [0, t]. To avoid numerical issues related to generation of observations too closely, the value of y * is set to t − t min if g(y) is negative for all t ∈ [0, t] or y * is too close to t, i.e., 0 < t − y * < t min .
After that, the value x k+1 is calculated from the value x k+1,① −1 substituting ① −1 by y * (without resolution of the ODEs and re-simulation of the system at t = t k+1 ). Once the value x k+1 is obtained, the observation t k+1 is computed starting from t k + ( t − y * ). Finally, the switch is performed by AC T I O N (varargs) and the simulation continues decreasing t = y * . More in detail, the function AC T I O N (varargs) is responsible for performing a discrete transition on the hybrid automata and, if required, updating the values of its discrete variables (e.g., the hybrid automata's state). Here, the function's parameter varargs holds all the information required to carry out the discrete transition (e.g., the current state q, the last computed variable x k , and the previously calculated variables x k−1 , x k−2 , ...,). At the next iteration, if zero-crossings are not detected then the default value t = t max is automatically restored (see Algorithm 1lines 14-20).
It is worth noting that the function g(y) can be negative for y ∈ [0, t] but no zero-crossings are present in such interval. This situation may occur when a zero-crossing was not correctly determined during the previous iterations for dif-ferent reasons. In this case, since by construction of g(x, q) the system works correctly only when g(x, q) ≥ 0, then the simulation can be incorrect at this interval and some additional actions are required. Since the value of t decreases with each iteration, the algorithm continues to simulate the hybrid systems with an ever smaller value of t up to the minimum allowed value of t min (see Algorithm 1-line 19), this behavior allows to stabilize the system after some iterations, but, in general, more actions can be required to solve this issue.

Hybrid systems used in the experiments
This section presents two well-known dynamic hybrid systems that exhibit the Zeno phenomenon, i.e., the Water tanks and Thermostat [see Johansson et al. (1999), Johnson et al. (2004)]. These systems have been considered in this work, since they are important both from practical and numerical points of view. More in detail, both the systems are linear hybrid systems, which means that the rate of change of each variable describing the system is constant, and the terms and conditions involved in the invariants, guards, and assignments are linear [see, e.g, Alur et al. (1995)].
Due to the Zeno phenomenon, the simulation of these systems fails because of numerical errors that start growing beyond the Zeno point. Such errors lead to incorrect simulation, and therefore, simulation results can be meaningless, especially when the systems are modeled with nonlinear functions.

Water tanks
Water tanks is an hybrid system that consists of two tanks t 1 and t 2 containing water [see Johansson et al. (1999)]. With reference to a generic Water tank t i with i ∈ I = {1, 2}, define x i as the water level, r i as the critical threshold, v i > 0 as the constant water flow going out. Define w as the constant flow of water that goes exclusively to either tank through a pipe at a given time t. Figure 2 shows the hybrid automaton model of the Water tanks system. The two blocks represent the two discrete states, whereas the physical behavior is described with the following differential equations: The objective of the system is to keep, for the tank t 1 the water level x 1 above r 1 , and for the tank t 2 the water level x 2 Fig. 2 Water tanks system hybrid automaton model above r 2 , assuming that x 1 (0) > r 1 ∧ x 2 (0) > r 2 . To pursue this objective while the water levels x 1 and x 2 keep dropping, when in one of the tanks, for example t 1 the water level drops below the critical threshold, i.e., x 1 < r 1 , the pipe switches to deliver the water, with the constant flow w, to the tank t 1 . As the water rapidly goes out of the tanks, the switching frequency of the pipe increases until it reaches the limit point that occurs when for each tank t i becomes , v 2 are positive constants, r 1 and r 2 are real constants in the following model: If w = v 1 + v 2 , then the system should stabilize on the level x 1 = r 1 and x 2 = r 2 , making infinitely many switches of the state starting from the time t c . If w < v 1 + v 2 , the system should also stabilize in one Water tank, but the level of the water in the other tank can decrease with time, making infinitely many switches as well. If w > v 1 + v 2 , then the system will not stabilize, and the water level in the tanks will grow. (In the latter case, the system is not Zeno, since the intervals between the switches will also grow.) The Water tanks hybrid system is also very important from a practical point of view. First, it has two different states q ∈ {0, 1} and for each fixed state q, there are two independent systems, i.e., Water tanks. Second, this system involves Zeno behavior already from the first occurrence g(x(t), q) = 0. In this case, the errors related to the Zeno phenomenon can arise very quickly, leading so to wrong states, incorrect simulation, and/or large oscillations of the water around the desired level r (r = r 1 and r 2 , respectively, for the first and the second tank). In the Water tanks system, the water level can be below the required level r . In this case, the simulation is correct only if the current state q leads to the growth of the water level in this Water tank. However, if the required water level r is not a constant, but a nonlinear function r = r (t), then there could be more than one Zeno point and the system's behavior becomes more sophisticated. In this case, it can be inefficient to generate observations too dense already from the first occurrence of g(x(t), q) = 0, and as will be shown below, the IC_RK4 method can be extremely efficient from the computational point of view.

Thermostat
A Thermostat is a digital component that senses the temperature of a room and performs actions so that the room's temperature x is always maintained near a desired set point σ by turning "off" and "on" a heater. When the Thermostat system starts, the heater is assumed to be "on" with an initial room temperature x, such that x < σ . In this situation, the room temperature increases according to the equationẋ = w−x, where w is the heater constant. The heating phase continues until x < σ. Upon the temperature reaches σ , the Thermostat turns the heater "off," and then, the room's temperature starts decreasing according to the equationẋ = −x.
Similar to the heating phase, the cooling one continues until x ≥ σ . When the temperature value becomes less than σ , the Thermostat turns the heater "on" and starts heating again. Note that, even if the Thermostat system looks simpler than the Water tanks one, it allows to better appreciate specific features of IC_RK4 method as the presence of the single threshold σ makes the handling of the Zeno phenomenon even more challenger [see Johansson et al. (1999), Johnson et al. (2004)]. Figure 3 shows the hybrid automaton model of the Thermostat system. The two blocks represent the two discrete states, i.e., "off" and "on." In each discrete state, the room's temperature x evolves according to the following differential equation: Fig. 3 Thermostat system hybrid automaton model Let x = x(t) ∈ R 2 , t ∈ [0, T ], q ∈ {0, 1}, w be a nonnegative real constant large enough, whereas r be a real constant, rand() generates a random number in the interval [0, 1]: Formula (10) describes the hybrid dynamics of a Thermostat for monitoring the internal temperature of a building. During the simulation, two outputs are provided, i.e., x and q that represent the temperature and the operation mode, respectively. In each mode, a specific differential equation regulates the temperature. The initial mode is "on" with an initial temperature value w, and the transition conditions between the discrete state q = 0 ("off") and q = 1 ("on") are regulated by AC T I O N (q).
The Thermostat system is also very important from the practical point of view, since it also has a clear physical realworld interpretation. Even if it is similar to the Water tanks system, it has several important differences with respect to it. First, it has only one variable x(t), instead of two independent systems x 1 (t) and x 2 (t), as in the Water tanks. Second, the function g(x, q) is randomized, so the zero-crossings are determined with a random error ξ(t), which can be considered as a standard white noise and can add some difficulties both for the proposed and standard methods. However, as will be shown below, the proposed algorithm can be useful in this case as well, showing promising behavior with respect to standard methods.

Description of Numerical experiments
All numerical experiments have been executed in MATLAB version R2016b. A software simulator of the Infinity Computer has been used for this purpose [see Falcone et al. (2020b) for details]. The initial and final times t 0 and T have been taken to 0 and 20, respectively, for all test problems. All parameters for each test problem have been set to the presented values only for simplicity: The obtained results and conclusions with other initial values and system parameters are similar to the presented ones.
The parameters of the Infinity Computer have been set for all test problems as follows. The simplest infinitesimal ① −1 has been used to set up the IC_RK4 method and the precision N has been set to 5, since the Runge-Kutta method of order 4 (RK4) method has the global error of order 4. For finding the zero-crossing points of g(y) for y ∈ [0, t], the standard bisection method up to machine precision has been used just for simplicity, since the variables x(t) are always monotonically decreasing or increasing around the zero-crossings g(x, q) = 0 in the studied test problems. (Other methods for the determination of zero-crossings based on global optimization algorithms from Casado et al. (2002), Molinaro and Sergeyev (2001) and Sergeyev et al. (1999) on more difficult real-life problems can be exploited.) Each of the hybrid systems presented above has the parameter(-s) r representing the simplest case for the zerocrossing function g: We want to keep the evaluations along the direction x = r (for the Water tanks system, there are two parameters r 1 and r 2 of the same meaning: x 1 = r 1 should be kept, if the current state is 1, and x 2 = r 2 , otherwise). However, in practice, the zero-crossing function g can have a nonlinear behavior: e.g., the cooling of the Thermostat along a nonlinear rule can be required. For this reason, for each test problem, we considered the following two series of experiments.
First, a simple case, when r is equal to a constant value ω, i.e., r = ω (r 1 = ω and r 2 = ω for Water tanks) is considered. Second, a more complicated case, when r = r (t) (r 1 = r 1 (t) and r 2 = r 2 (t) for Water tanks), is considered, where r (t), r 1 (t), r 2 (t) are nonlinear functions.
The already mentioned method RK4 has been used for solving ODEs. For each test problem, first, the standard RK4 method with the fixed stepsize t has been applied (we will call this method as RK4_F, hereinafter). Then, the IC_RK4 method described above with the dynamic stepsizes and parameters t min and t max is applied for the same problems.
It should be noticed that in this paper, the internal MAT-LAB methods with dynamic stepsizes (e.g., Dormand-Prince methods) are not considered for several reasons. First, these methods are just related to different numerical algorithms for solving ODEs and not to the methods of generating observations t k : These methods do not generate the observations dynamically, but use more sophisticated numerical methods for solving ODEs, with respect to the standard RK4 method, thus, comparing them with the RK4 method is not correct.
Moreover, these methods are not adapted well for the case, when g(x, q) is nonlinear, allowing only to determine simple zero-crossings of a type x = 0. However, it should be also noted that the algorithm for generating the observations proposed in this paper can be used with the above-mentioned dynamic methods, as well.
For each test problem, the standard method RK4_F with the fixed stepsize t has been applied for simulating the system using different values of t: 0.5, 0.05, 0.005, 0.0005, and 0.00005. Then, the IC_RK4 method has been applied for the same problems using three different values of the parameters t min and t max : [ t min , t max ] = [0.005, 0.5], [0.0005, 0.05], and [0.00005, 0.005] in order to compare the obtained results with the standard method RK4_F with the stepsizes t of the same ranges. In the following subsections, only the figures related to the stepsizes t = 0.5, 0.05, and 0.005 for RK4_F, and [ t min , t max ] = [0.5, 0.005] for IC_RK4 are presented; to improve the readability of the paper, the figures with the simulation results of other stepsizes are presented in Appendix (avoiding a lot of additional figures in the main text). This is done because these additional figures show the same behavior of the algorithms, but they substantiate the conclusions presented in this section, so they should be also attached to the paper. The numbers of generated observations by both methods are presented in each figure. They are also summarized in Tables 1, 2 for the sake of completeness.

Case r 1 = r 2 = !
In this section, the test problem studied numerically is the Water tanks system described in Sect. 3.1. The following parameters have been set for this system: x(t 0 ) = [x 1 (t 0 ), x 2 (t 0 )] = [2, 0], the initial state q = 0, and the parameters w, v 1 , and v 2 have been chosen in the way that w = v 1 + v 2 (since it is required to keep the water levels x 1 and x 2 at the constant levels r 1 and r 2 : w = 5, v 1 = 2, v 2 = 3). The desirable water levels r 1 and r 2 for both the Specifically, for the Water tanks system with r 1 = sin(t) − t and r 2 = sin(t), and for the Thermostat system r (t) = ω + sin(t) − t with ω = 70 Fig. 4 Simulation of the Water tanks system performed with the standard method RK4_F with ω = 1; thus, r 1 = r 2 = ω and the fixed stepsizes t = 0.5. Simulation results show that the standard method can lead to wrong states of the systems and its incorrect behavior, which can be resolved only using smaller t. Generated observations t k , k = 0, ..., N , with N = 40) are dense everywhere (the observations are indicated below the graphs by signs "+") Water tanks have been set to a constant value ω = 1, i.e., r 1 = r 2 = 1. Results of the simulation by RK4_F using t = 0.5, 0.05, and 0.005 are presented in Figs. 4, 5, and 6, results of Fig. 5 Simulation of the Water tanks system performed with the standard method RK4_F with ω = 1; thus, r 1 = r 2 = ω and the fixed stepsizes t = 0.05. Also in this case, the standard method RK4_F leads to incorrect system's behavior. Generated observations Fig. 6 Simulation of the Water tanks system performed with the standard method RK4_F with ω = 1; thus, r 1 = r 2 = ω and the fixed stepsizes t = 0.005. The same considerations, which were described in both the cases t = 0.5 and t = 0.05 about the results of the standard method, can be done here. Generated observations t k , k = 0, ..., N , with N = 4000 are dense everywhere (the observations are indicated by signs "+") Figures 4, 5, and 6 show that the standard fixed stepsize method RK4_F also can lead to an incorrect system's state and, as a consequence, to incorrect behavior. In particular, Fig. 4 shows the incorrect behavior of the Water tanks system in the interval [0, 0.5]. More in detail, in Tank 1 the water level decreases and remains below the desired value r 1 , but it should increase to reach and stay on r 1 (see Fig. 4a). On the other tank, i.e., Tank 2, the water level increases above the desired value r 2 , while it should decrease to reach and remain on r 2 . However, with smaller values t, the method converges to the correct simulation (see Figs. 5,6).
In its turn, the IC_RK4 method (see Fig. 7) simulates the system similarly to the standard method with the smallest value t = 0.005, switching right away the system from the "wrong" initial state to the correct one at the first iteration. Again, since the IC_RK4 method generates the observations Fig. 7 Simulation of the Water tanks system performed with the evaluated method IC_RK4 with ω = 1; thus, r 1 = r 2 = ω and t min = 0.005, t max = 0.5. The simulation results show that the proposed IC_RK4 method with the dynamic stepsizes t and automatic zero-crossing checking allows to keep the same precision as the fixed stepsize method with small t, but generates less observations. Generated observations t k , k = 0, 1, ..., 3188, are dense only after the first zero-crossing (the observations are indicated below the graphs by signs "+") t k more frequently, where it is required, and less frequently on "good" and "stable" regions, then it generates fewer observations obtaining at least the same accuracy that the standard fixed stepsize method does. One can see also that the IC_RK4 method generates always fewer evaluations keeping at least no worse precision of the simulation. Both Table 1 and Fig. 7 show that the IC_RK4 method starts to generate the dense observations only after the first zero-crossing. This allows to use fewer computational resources on the initial subintervals where there are no zero-crossings and the system is well conditioned.
However, it should be noted that the IC_RK4 method does not resolve completely the Zeno phenomenon: It also can lead to incorrect behavior beyond the Zeno points. In Fig. 7,  Fig. 8 Simulation of the Water tanks system performed with the standard method RK4_F with r 1 (t) = sin(t) − t, r 2 (t) = sin(t), and fixed stepsizes t = 0.5. Simulation results show that the standard method leads to wrong states and incorrect system's behavior, which can be resolved only using smaller t. Generated observations for example, the water level in the second Water tank can "jump" a little below the level r 2 at some moments t k (at the same moments, the water level in the first tank "jumps" above the level r 1 ), this is related only to the Zeno behavior of the system and is due to Algorithm 1, lines 14-15: At some moments, the zero-crossing point y * becomes too close to t, so it is replaced by t − t min to avoid numerical issues related to ill-conditioning, which leads to these "jumps."

Case r 1 (t) = sin(t) − t and r 2 (t) = sin(t)
A numerical experiment has been conducted by changing the parameters r 1 and r 2 to nonlinear functions r 1 (t) = sin(t) − t and r 2 (t) = sin(t). The functions r 1 (t) and r 2 (t) have been chosen to study the case when the Water tank as should behave differently: We want to decrease the water level x 1 Fig. 9 Simulation of the Water tanks system performed with the standard method RK4_F with r 1 (t) = sin(t) − t, r 2 (t) = sin(t), and fixed stepsizes t = 0.05. Also in this case, the standard method leads to incorrect simulation of the system. Generated observations t k , k ∈ [0, ..., N ], with N = 400 are dense everywhere (the observations are indicated with signs "+") in the first tank along with the level r 1 (t), while we would like to keep the level x 2 in the second tank higher than r 2 (t). For this reason, the velocities of the flows w, v 1 , and v 2 have been fixed in the way that w < v 1 + v 2 .
Also in this case, the standard method RK4_F was able to simulate the system behavior only using smallest values of t (see Figs. 8,9,and 10), while the proposed one successfully simulated the system by generating much fewer observations than the standard one (see Fig. 11). One can see that the IC_RK4 method has generated more dense observations only in "unstable" and ill-conditioned regions. Table 2 reports the number of observations performed by the two methods to Fig. 10 Simulation of the Water tanks system performed with the standard method RK4_F with r 1 (t) = sin(t) − t, r 2 (t) = sin(t), and fixed stepsizes t = 0.005. Same considerations, which were described in both the cases t = 0.5 and t = 0.05 about the adoption of the standard method to simulate the system. Generated observations t k , k ∈ [0, ..., N ], with N = 4000 are dense everywhere (the observations are indicated with signs "+") simulate the system with r 1 (t) = sin(t) − t and r 2 (t) = sin (t).
Also in this case, the advantages of the IC_RK4 method become more evident in the case of a nonlinear function g(x, q): It has generated almost 14, 34, and 62 times less observations than the standard one for t = t min = 0.005, 0.0005, and 0.00005, respectively (282 vs 4001, 1156 vs 40,001, and 6387 vs 400,001 observations, respectively, see Table 2).

Case r = !
The example described in Sect. 3.2 is different with respect to the previous one due to the presence of random perturba- Fig. 11 Simulation of the Water tanks system performed with the evaluated method IC_RK4 with r 1 (t) = sin(t) − t, r 2 (t) = sin(t), and t min = 0.005, t max = 0.5. The simulation results show that the IC_RK4 method with the dynamic stepsizes t and automatic zerocrossing checking allows to maintain the correct behavior of the system. Generated observations t k , k = 0, 1, ..., 281, are more dense around zero-crossings (the observations are indicated below the graphs by signs "+") tions. It has been simulated using the following parameters: x(t 0 ) = x 0 = 60, initial state q = 0, w = 100, the desirable temperature r has been set to a constant value ω = 70 o F, random numbers have been generated on the interval [0, 1] using the internal MATLAB function rand() (the seed of the generator has been set to 1 by the command rng(1)).
Results of the simulation by RK4_F using t = 0.5, 0.05, and 0.005 are presented in Fig. 12, results of IC_RK4 using [ t min , t max ] = [0.5, 0.005] are presented in Fig. 13. The remaining results are presented in Appendix (see Figs. 20,21). Figures 12 and 13 show that the standard method RK4_F leads to an incorrect simulation of the system with biased mean temperature: it is required to keep the temperature of the Thermostat at the level of 70 o F, while one can see that The simulation results show that the standard method RK4_F can lead to wide oscillations of the system's variables. These oscillations can be stabilized only using smaller t, and, therefore, can lead to biased average temperature with respect to the desirable level r o F. Generated observations t k , k = 0, ..., N , (N = 40, 400, and 4000, respectively for a-c) are dense everywhere (the observations are indicated below the graphs by signs "+") the mean level of the temperature generated by RK4_F is below this level (this is avoided only using smaller t, see Unlike the standard method RK4_F, the IC_RK4 one with the dynamic stepsizes t and automatic zerocrossing checking allows improving the accuracy of the simulation decreasing the oscillations and keeping them around the desirable level r o F. Generated observations t k , k = 0, ..., 850, are more dense after the first zero-crossing (the observations are indicated below the graphs by signs "+") the respective figures in Appendix). The simulation behavior in Fig. 20 seems to be deterministic just because the value t, in this case, is too large and therefore dominates the stochastic perturbations. On the other hand, the IC_RK4 method simulates the system correctly generating the temperature of the desirable level on average (see Fig. 13). Table 1 shows that the proposed algorithm has generated significantly less observations than the standard method RK4_F: almost 5, 31, and 77 times less for t = t min = 0.005, 0.0005, and 0.00005, respectively (851 vs 4001, 1266 vs 40,001, and 5157 vs 400,001 observations, respectively, see Table 1).
It should be noted that the IC_RK4 method determines the zero-crossings in this problem with a random error, since the function g(x, q) contains random numbers. Thus, oscillations of the variable x(t) with smaller t max and t min have higher amplitudes with respect to the standard method RK4_F (see Fig. 21). However, first, as for the standard method, smaller values of t min and t max allow reducing these amplitudes. Second, the proposed algorithm is able to maintain a correct simulation in this case as well as for the previous hybrid systems, decreasing the number of observations significantly: e.g., as was mentioned above it has generated almost 77 times fewer observations using [ t min , t max ] = [0.00005, 0.005] than the standard method using t = t min = 0.00005.

Case r = ! + sin(t) − t
A further numerical experiment has been conducted by changing the parameter r = ω to r (t) = ω + sin(t) − t with ω = 70 similarly to the previous experiments. From the Fig. 14 Simulation of the Thermostat system performed with the standard method RK4_F with the constant value ω = 70; thus, r (t) = ω + sin(t) − t and t = 0.5 (a), t = 0.05 (b), t = 0.005 (c). The simulation results show that the standard method RK4_F can lead to wide oscillations of the system's variables, is able to stabilize them only using smaller values of t, and can lead to biased cooling. Generated observations t k , k ∈ [0, ..., N ], with N = 40, 400, and 4000, respectively for a-c) are dense everywhere (the observations are indicated below the graphs by signs "+")

Fig. 15
Simulation of the Thermostat system performed with the evaluated method IC_RK4 with the constant value ω = 70; thus, r (t) = ω + sin(t) − t and t min = 0.005, t max = 0.5. Unlike the standard method RK4_F, the IC_RK4 method with the dynamic stepsizes t and automatic zero-crossing checking allows to improve the accuracy of the simulation decreasing the oscillations and keeping the desirable cooling. Generated observations t k , k = 0, 1, ..., 977, are dense after the first zero-crossing (the observations are indicated below the graphs by signs "+") physical point of view, this case can be described as follows. At the time t = 0, the required temperature is ω o F. After that, a cooling along with a nonlinear rule r (t) is required (for example, a too fast or slow cooling can damage the simulated devices).
The results of the simulation by RK4_F using t = 0.5, 0.05, and 0.005 are presented in Fig. 14, the results of IC_RK4 using [ t min , t max ] = [0.5, 0.005] are presented in Fig. 15. The remaining results are presented in Appendix (see Figs. 22,23).
In this case, the simulation results confirm that the standard method was able to simulate the system behavior correctly only using the smallest values of t (see Fig. 14), while the IC_RK4 method successfully simulated it by generating significantly fewer observations (see Fig. 15).
The obtained results confirm the advantages of the IC_RK4 method: it is able to generate fewer observations with respect to the standard one simulating correctly the desirable cooling. In particular, the IC_RK4 method has generated almost 4, 29, and 75 times less observations for t = t min = 0.005, 0.0005, and 0.00005, respectively (978 vs 4001, 1403 vs 40,001, and 5314 vs 400,001 observations, respectively, see Table 2). However, again, for the same reasons as in the previous case, i.e., due to random errors in the determination of zero-crossings, the IC_RK4 method simulates the system with a higher amplitude with small t min and t max , than the standard method, maintaining always a correct behavior of the system and generating a significantly smaller number of observations with respect to the standard method.

Conclusion
To support the design, development, and operation of modern Cyber-Physical Systems, many research efforts are focusing on the definition of methods, models, and techniques capable of capturing the interactions between the physical and cyber components through the definition of hybrid system models. Unfortunately, one of the most important issues in hybrid systems is the Zeno phenomenon consisting in the identification of zero-crossings represents a crucial aspect to adequately simulate these systems. The numerical methods adopted by classical modeling and simulation techniques are slow because they require solving the ordinary differential equations at each step. To overcome this issue, and, therefore, avoid redundant computations, the IC_RK4 method working numerically with finite, infinite, and infinitesimal on the Infinity Computer has been used.
The IC_RK4 method has been exploited to study two well-known Zeno hybrid systems, i.e., Water tanks and Thermostat. Both the systems have been modeled with constant and nonlinear threshold functions to stress the IC_RK4 method, and evaluate its performance in detecting zerocrossings by generating time observations using the infinite quantity ① offered by the Infinity Computer, dynamically.
Simulation results have confirmed the validity and performance of the IC_RK4 method in simulating Zeno hybrid systems. Moreover, it has been shown in the performed experiments that the IC_RK4 method allows improving the accuracy of the simulation, also in the presence of nonlinear threshold functions, because the effective detection of zerocrossings directly influences the efficiency of the simulation.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.
Data Availability All data supporting the results reported in the article are present in the paper.

Declarations
Conflict of interest All authors declare that they have no conflict of interest.
Human and animal rights. This article does not contain any studies with human participants or animals performed by any of the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Water tanks system
See Appendix Figs. 16, 17, 18, and 19.   Fig. 16 Simulation of the Water tanks system performed with the standard method RK 4 F with the constant value ω = 1; thus, r 1 = r 2 = ω and the fixed stepsizes t = 0.00005. Trend of the water level related to Tank 1(a), Tank 2(b) Fig. 17 Simulation of the Water tanks system performed with the evaluated method I C RK 4 with the constant value ω = 1; thus, r 1 = r 2 = ω and [ tmin, tmax] = [0.0005, 0.05]. Trend of the water level related to Tank 1(a), Tank 2(b) Fig. 18 Simulation of the Water tanks system performed with the standard method RK 4 F with r 1 = sin(t) − t, r 2 = sin(t) and fixed stepsizes t = 0.00005. Trend of the water level related to Tank 1(a), Tank 2(b) Fig. 19 Simulation of the Water tanks system performed with the evaluated method I C RK 4 with r 1 = sin(t) − t and r 2 = sin(t) and