# Qualitative tendencies for hybrid system simulation

Baptiste Gueuziec, Frédéric Boulanger, Jean-Pierre Gallois

► To cite this version:

Baptiste Gueuziec, Frédéric Boulanger, Jean-Pierre Gallois. Qualitative tendencies for hybrid system simulation. MPM4CPS 2023 - Multi-Paradigm Modelling for Cyber-Physical Systems -26th International conference on model driven engineering languages and systems, Oct 2023, Västerås, Sweden. hal-04188249

HAL Id: hal-04188249

https://hal.science/hal-04188249v1

Submitted on 25 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Qualitative Tendencies for Hybrid System Simulation

Baptiste Gueuziec

Université Paris-Saclay

CEA, List

F-91120 Palaiseau, France

baptiste.gueuziec@cea.fr

Jean-Pierre Gallois

Université Paris-Saclay

CEA, List

F-91120 Palaiseau, France

jean-pierre.gallois@cea.fr

Frédéric Boulanger

Université Paris-Saclay, CNRS, CentraleSupélec

Laboratoire Méthodes Formelles

F-91190, Gif-sur-Yvette, France

frederic.boulanger@centralesupelec.fr

**Abstract**—The study of complex systems is essential in many scientific fields for safety, reliability, and security reasons. These systems are increasingly complex and often hybrid (with continuous behaviors in discrete modes), making modeling tasks more difficult. Each design requires a tradeoff between acceptable computation time and good precision. Simulations using numerical methods often sacrifice time efficiency for more precise results, but this choice can show its limits on complex systems. Qualitative reasoning techniques were proposed to sacrifice precision and knowledge to improve speed and generality. They are especially used in diagnosis and fault isolation. Eventually, abstraction methods such as interval propagation or flow-pipe execution emerged in what we now call semi-qualitative simulation, mainly used for precision and property proving. This article proposes an approach that bridges qualitative analysis and numerical simulation methods. It consists in building a qualitative model of a system that is used as a map of its behavior in various regions of its state space to improve the efficiency of its simulation. Therefore, the joint use of models of different natures helps analyze the system's behavior.

**Index Terms**—Cyber-physical systems, simulation, qualitative analysis

## I. INTRODUCTION

Complex systems [1] are involved in the modeling and analysis of multi-component systems and natural phenomena. Running simulations of such systems requires an efficient representation of the model and a robust methodology. Modeling those systems implies isolating their most important characteristics and optimizing their representation for computation. A significant difficulty is to find the best complexity/precision trade, especially for hybrid systems [2], which need to consider both discrete and continuous behaviors.

Depending on the available knowledge and the goal of the simulation, the modeling phase will use different methods to abstract the system and represent the information. Classically, more precise representations using numerical data and differential equations are used for continuous behaviors. The simulation consists in the numerical integration of these equations using methods such as Euler or Runge-Kutta to compute the precise behavior of the system.

However, Brown [3] and De Kleer [4] introduced the idea of a qualitative representation of knowledge using symbolic expressions and relations between the components of a model. They developed it to help in computer-assisted calculus and

electronic problem-solving. Since then, scientists such as Hayes [5], [6] and Forbus [7] opened new perspectives by introducing notions such as naive physics and adding the notion of process and dynamics. The question of high-scale knowledge representation of a complex environment leads to ideas such as concept-based graphs. Meanwhile, the challenge of dynamics representation permitted the use of causality graphs [8] or bond graphs [9]. Kuipers studied the application of qualitative modeling for simulation [10], [11] and proposed tools to simulate models with qualitative methods (QSIM) using qualitative differential equations or semi-qualitative techniques (Q3) using interval propagation. He represents models using qualitative differential equations, which are abstractions of ordinary differential equations (ODE) obtained by isolating the extrema, the variation sense, and the qualitative value of the variables. The evolution of the values is represented by qualitative states separated by landmarks (local extrema of one variable or one of its derivatives). His tool Q3 combined this approach with interval propagation techniques.

More advanced simulation methods such as the works on flow-pipe propagation [12] or closed expansivity study [13] tried to correct some drawbacks of the propagation of intervals. Finally, Tiwari opened another branch of qualitative modeling in [14], more adapted to systems represented by ODEs. He proposed a numerical analysis of a system's polynomial equations and their derivatives and used the results to create a map of the state space of the model. Works presented in [15]–[17] studied these different approaches of qualitative reasoning and knowledge propagation, seeking an optimal representation of models for general tasks. In this paper, we propose a modeling method to improve simulation performance where precision is not the main criterion, and the goal is to identify qualitative behaviors without needing quantitative precision. We will use the different concepts of qualitative modeling, enriched with *qualitative tendencies*, which is the contribution that we will detail further.

## II. QUALITATIVE MODELING AND ORDERS OF MAGNITUDE

### A. *Dynamic System Formalism*

Throughout the rest of this document, we will consider the study of dynamical systems expressed as $D = \langle Q, X, F, I, T \rangle$,
<figure><img src="image_2.png" /><figcaption>Fig. 1. Hybrid model of a thermostat</figcaption></figure>

with $Q$ the set of discrete variables of the system. Each variable in $Q$ represents the mode of a part of the system (for instance, in a car, a variable may give the state of the cruise control, and another the state of the air conditioning). These variables take value in $Q$, which is the cartesian product of the domains of each variable. $X$ is the set of continuous variables, which have values in $\mathbf{X}$. $\mathbf{X}$ is classically $\mathbb{R}^n$, with $n$ the number of continuous variables. $F$ gives the flow of the system (a function associating each mode to a set of ODE on $X$ and its derivatives constraining the evolution of the continuous variables). $I$ is the function associating each mode to its invariant conditions. These invariants are properties that always hold in a mode, and that are useful to analyze the possible behaviors of the system. We represent $F$ and $I$ with dictionaries. $T$ represents the discrete transitions between modes as a set of quadruplets $(m_1, \text{cond}, m_2, \text{reset})$, which means that when the mode of the system is $m_1$ and condition $\text{cond}$ on the continuous variables is satisfied, the system goes into mode $m_2$ and resets the values of the $X_i$ using the function reset. The validity of $\text{cond}$ forces the transition.

For example, a thermostat can be represented as $Th = \langle Q, X, F, I, T \rangle$ [2], [14] with $Q = \{\text{off}, \text{on}\}$, $X = \mathbb{R}$, and:

$$Q = \{\text{mode}\} \text{ with } \text{mode}_{\text{init}} = \text{off},$$

$$X = \{T\} \text{ with } T_{\text{init}} = 75,$$

$$F = (\text{mode} = \text{on}: \dot{T} = -T + 100, \text{mode} = \text{off}: \dot{T} = -T),$$

$$I = (\text{mode} = \text{off}: T > 50, \text{mode} = \text{on}: T < 100),$$

$$T = \{(\text{on } \xrightarrow{T \ge 80} \text{off}, T \mapsto T), (\text{off } \xrightarrow{T \le 70} \text{on}, T \mapsto T)\}$$

It can be represented as an hybrid automaton (see Fig. 1).

When speaking about the equations of a model, we will include every expression in the flow, the invariants, and the conditions of the transitions. We will limit our study to polynomial (including multivariate polynomials) functions. We are currently working on generalizing our results to logarithmic-exponential functions (a family of functions closed under sum, product, and composition) and to more general functions.

For the qualitative analysis of the system, we need to introduce the matrix $q$ of the characteristic values of the system. This matrix of size $|Q| \times |X|$ gives the characteristic value of each continuous variable for each mode of the system. The characteristic value $q_{i,j}$ gives us information about the nominal range of the variable $X_j$ in the mode corresponding to the $i^{th}$ element of $Q$. In our example, we could define $q = \begin{pmatrix} 75 \\ 75 \end{pmatrix}$ considering off as mode 0 and on as mode 1. We

<figure><img src="image_3.png" /><figcaption>Fig. 2. Qualitative model of the thermostat</figcaption></figure>

chose these values because $T$ varies in both modes mainly in the interval $[70, 80]$.

## B. Qualitative Modeling of Dynamic Systems

The dynamics of complex systems is usually described with modes, states, and events [7]. A mode corresponds to the vector formed by the value of each discrete variable in $Q$ and a change of mode can modify the flow [14] and the causal ordering of the events in the system [8]. The state of the system is described using both the mode and the numerical value of each continuous variable in $X$, and corresponds to a precise situation. The drawback of using numerical values is the lack of generality: simulating the system with specific numerical values for continuous variables requires either studying one trajectory or using many computational resources to allow interval propagation. However, it tells us nothing about the general qualitative behavior of the modeled system. For this reason, in qualitative reasoning, we choose to avoid the computation and propagation of precise values. Instead, we focus on the relative position of both variables and their derivatives compared to significant thresholds or frontiers defined by equations. These equations allow us to introduce the concept of qualitative states. The relative position of the current state of the system to these frontiers defines the qualitative state of the system at an instant. An example of this qualitative division is shown in Fig. 2, which is a qualitative representation of our example thermostat model. The black rectangles correspond to the modes of the system (here on and off). The black horizontal arrows correspond to the transitions between modes. The circles labeled in blue correspond to qualitative states (increasing, decreasing, min, and max), and the blue vertical transitions between these states
are triggered according to qualitative thresholds computed
from the equations of the system.

### C. Orders of Magnitude

Among the various tools used to reason on qualitative
models, we can find the orders of magnitude [18]. This method
expands the reasoning on sign, which is the very base of
qualitative methods [7], [10]. It is mainly supposed to avoid
the limitations of the sign algebra, especially the difficulty of
applying the + (and less significantly ×) operators to values
of unknown magnitude. Taking the magnitude of a value into
account allows for more precise and coherent computations but
also has drawbacks on the permissiveness of operators [19].
For example, transitivity, associativity, and distributive prop-
erties are more limited in this qualitative algebra. Reasoning
about magnitude can be absolute (comparison to fixed values)
or relative (comparison between variables) [20].

The first option to reason about orders of magnitude implies
partitioning the state space according to constant landmarks.
Variables are individually compared to these reference val-
ues $k_i$ and placed in the corresponding rectangle (multi-
dimensional closed interval), determining their magnitude. It
also requires choosing a scale to apply to these reference
values that will set the abstraction's regularity and granularity.
The interest is that by optimizing the choice of landmarks,
it is possible to abstract the precise values completely while
preserving helpful information about the value of the state
variables. The scale can, for example, be linear or, if the
variable is to stay positive, logarithmic and rely on powers
of $k_i$. The values $k_i$ will here depend on design criteria
of the system, such as the initial value or the anticipated
extreme points. When the sign can change, a linear scale may
be preferred. More generally, the choice of the state space
partition will strongly depend on the simulation's objectives
and the system's design [21]. For example, a linear scale in one
dimension based on a unique value $k$ would be structured with
landmarks on the values 0, $-k$, $2k$, $-2k$, ... It could be used
in the case of the model of the thermostat, where the unique
system variable varies between values that are simultaneously
closely spaced and sufficiently far from zero to make any
logarithmic scale unusable. In this case, an adapted scale could
be a linear one with a value $k = 2$.

A second option is to use relative orders of magnitude. This
option compares variables with each other rather than with
reference values. This requires the physical quantities and units
to be comparable. Different comparison systems exist, such
as FOG [20], O(M) [18], or Rom [22]. These systems do not
use the same operators, but some of them are used in most
systems: Ne (negligible, also noted ≪), Co (comparable),
and Vo (Neighbor, or very close to). These operators can
give more precision than the usual ones, such as =, <, >, or
≈. Orders of magnitude are mainly used in purely qualitative
algebra because the intrinsic properties are incompatible with
the archimedean property and therefore not supported in $\mathbb{R}$.
For example, it is supposed in Rom that if a Ne b, therefore
$(a+a)$ Ne b, which is not correct in archimedean spaces. Some

recent works added the operator $Di$, meaning distant [23],
and corresponding to the situation where the ratio between
two values is too important to consider them as comparable
(not of the same order of magnitude), but too small for
one of them to be negligible. Absolute and relative orders
can also be combined, but this requires an adapted space
partition and a precise definition of the relative operators.
The difficulty is to make the different frontiers match each
other. Otherwise, we get the superposition of two independent
measure scales, adding too much complexity for few benefits.
These techniques were designed to compare values, but we
must also compare functions to compare the behaviors of
variables and not only static values.

### D. Order of Growth

Many works of the previous century have studied categories
of functions, and some categorized them on their order of
growth. This concept eases the representation of the behavior
of a function $f : t \mapsto f(t)$ when $t \to \infty$. The notion of
order of growth was first formalized by Borel [24] and pushed
further by other scientists using the works of Hardy [25] on
logarithmic-exponential functions (set of functions obtained
by addition, multiplication, or composition of logarithm or
exponential functions, closed under addition, multiplication,
and composition). Let us consider $f : \mathbb{R} \mapsto \mathbb{R}$ a function. We
define the order of growth of $f$ by $c(f) = \lim_{t \to +\infty} \frac{|\ln|f(t)|}{\ln(t)}$ and
$c(f : t \mapsto 0) = -\infty$. In the case of polynomial
functions $f : t \mapsto \sum_{i=0}^{n} a_i t^i$, we have $c(f) = n$. When $f$
is the sum of many terms, $c(f)$ only considers its asymptotic
dominant term without consideration of the others. Moreover,
$c(f)$ only has a meaning if $f$ converges or diverges to an
infinite: if $f$ is periodic or has no clear limit when $t \to \infty$,
then $c(f)$ does not exist. The expression of $c$ allows us to
highlight some properties of the order of growth regarding its
reaction to operations. Using the properties of the logarithm, it
appears that for two logarithmic-exponential functions $f$ and
$g$, $c(fg) = c(f) + c(g)$ and $c(f \circ g) = c(f)c(g)$. Moreover,
if $c(f) \neq c(g)$, then $c(f + g) = \max(c(f), c(g))$. Still, this
measure of asymptotic behaviors has some weaknesses. First,
with the given formula, it is impossible to distinguish constant
functions from logarithmic ones or to make a difference
between two exponential functions with their order of growth.
Actually, $c(f \mapsto k \neq 0) = c(f : t \mapsto \ln(t)) = 0$, and also
$c(f : t \mapsto e^t) = c(f : t \mapsto 10^{2t}) = \infty$, which is a severe
limitation for a behavior classification. Classifying an infinite-
diverging and a constant function with the same value cannot
be sufficient as a characterization. Finally, for our usage, this
order of growth only applies when $t \to \infty$, making it useless
for analyzing systems at finite time scales through simulation.
Therefore, we want to generalize this concept to local areas
of the state space.

### III. QUALITATIVE TENDENCIES

Our objective is to extend the notion of order of growth
to local rectangles to offer a qualitative abstraction of the
functions of a system using piecewise decomposition. To this
end, we will combine the tools presented earlier to propose a qualitative local categorization and simplification of the equations of a model. We will then reduce the complexity and make the model's behavior easier to understand at the price of a loss in precision.

Let us consider the evolution of a variable $y$ that varies according to the function $y = (t-3)^3 - 4 * t^2 + 2 * t$. The observed behavior of the variable will be that of a classic cubic function, with two phases of polynomial growth separated by a decay phase. On a small scale, this decay period is essential to determine a stop in the increasing qualitative behavior observed earlier. However, on a larger scale, if we are just interested to know the difference of magnitude at two times $t_1$ and $t_2$, this decrease period has no interest: the tendency is that the variable has followed a cubic growth between $t_1$ and $t_2$. This difference between general behavior and local trajectory is the center of our idea and work. Depending on the level of granularity and the information we seek, we can adapt the expression of the equations of a system to optimize its study without altering the results.

## A. Comparison Tools and Dominant Terms

To generalize the notion of order of growth to a local interval, we also use the notion of order of magnitude (both absolute and relative) to simplify the behavioral equations of the model. The logic is, for a vector of variables $X$ and $f: X \mapsto f(X) \in \mathbb{R}$, to extend orders of magnitude from $f(X)$ to $f$ itself. By extension, this would apply to differential equations as well. Intuitively, a function $f(x) = x^2 + x + 1$ will not have the same kind of behavior if $x = 0.025$ or if $x = 1500$. In one case, the behavior of $f$ will (locally) be linear, while in the second, it will be polynomial with a major predominance of the $x^2$ element. More specifically, when $0 < x \ll 1$ (using the comparison operators introduced in the previous parts), it appears that $f(x) \approx 1$ because $x$ and $x^2$ are negligible compared to 1 in that area. Introducing the operator $c_{local}(f)$, we would have $c_{[-\epsilon, +\epsilon]}(f) = 0$, with $0 < \epsilon \ll 1$. As $\lim_{x \to 0} f(x)_{x \to 0} = \lim(1)_{x \to 0}$, it is possible to neglect the other terms of the function $f$ to obtain an acceptable approximation in the neighborhood of 1.

If we can give a more situational and larger negligibility criterion (for example, if we consider that in the studied system, $x \ll x_s$ if $x/x_s < k$, with $k$ a constant to define and $x_s$ a specific value for the variable $x$), then we could replace the neighborhood of finite values by a rectangle inside which $f(x) \approx 1$. On the opposite, we could also say that if $x > \frac{1}{k}$, then $f \approx x^2$ for the same reason. Of course, the value of $k$ must be chosen accordingly to the design of the system, and there is no question of giving a general value that could be applied anywhere. An idea to correctly choose $k$ is to fix the limits of $x$ (either explicitly with an invariant of the system or implicitly by reasoning on the nature and quantities of the variable) and to deduce a negligibility criterion from these limits. For example, a chemical system, with density variables evolving between 0 and 1, will not have the same negligibility criteria as in a solar-system model, with distances outgrowing

<figure><img src="image_4.png" /><figcaption>Fig. 3. Approximation of $f(x) = x^3 - 2x^2 - 2$ by dominance areas.</figcaption></figure>

billions of kilometers. Knowing the expected variation space of the simulation, we can deduce the highest magnitude that, if added to a value in this interval, would not change its interpretation and seem negligible compared to it. It will also depend on the allowed margin of error: the authorized inaccuracy will not be the same for a soap factory and a nuclear plant. An adaptation to every case may be necessary.

When $x \to \infty$, we use what we explained earlier regarding the growth of the sum of functions. As the order of growth is absorbent with the sum of monomials of different degrees, we consider that for each monomial $m$ of $f$ such that $c(m) < c(f)$, there exists a specific value $x_s$ (depending on the negligibility criteria mentioned earlier), such that for all $x > x_s$, $m(x) \ll f(x)$. We name $f_p$ a target function aimed to approximate $f$ precisely enough to preserve the properties of $f$ but with simpler behavior and expression. $f_p$ can be created by replacing locally $f$ by approximations in specific rectangles from specific landmarks by suppressing one after the other the monomials of $f$ when they become negligible compared to the dominant term (the highest degree monomial) of the function. We give an example of a simple function in Fig. 3, showing an approximation of the function $f: x \mapsto x^3 - 2x^2 - 2$ to a piecewise defined function.

In this figure, the two curves are mingled (excepted in the $|x| > 10$ areas) because our approximation $f_p$ (black curve) using the dominant terms in the different areas of the state space gives us a satisfying representation of the function $f$ (the red curve) in terms of approximation and evolution form. $f_p$ is defined as (the value 10 appears because we consider orders of magnitude):
$$f_p = \begin{cases} -2 & \text{if } |x| \le \frac{1}{\sqrt{10}} & x^3 \text{ and } x^2 Ne 2 \\ -2x^2 - 2 & \text{if } \frac{1}{\sqrt{10}} \le |x| \le \frac{1}{\sqrt{10}} & x^3 Ne 2 \text{ and } x^2 \\ x^3 - 2x^2 - 2 & \text{if } \frac{1}{\sqrt{10}} \le |x| \le \sqrt{10} & x^3 - 2x^2 \text{ if } \sqrt{10} \le |x| \le 10 \\ x^3 - 2x^2 & \text{if } \sqrt{10} \le |x| \le 10 & 2 Ne x^3 \text{ and } x^2 \\ x^3 & \text{if } |x| > 10 & x^2 \text{ and } 2 Ne x^3 \end{cases}$$

## B. State Space Sub-Partitioning

The point is now to apply this partition based on order of growth to study the different equations of the model and to map their behavior in the state space. As mentioned earlier, the approximated behavior of a variable depends on the scale but also on the time of observation. To know which abstraction will be authorized, we must first confront the equations to the magnitude taken by the variable in the state space and then study the time slot we are interested in. This will require a mapping of the state space that associates to each set of the partition the dominant terms to be kept in the approximated equation. To take our first example back, the dominant term of the function $x \mapsto x^2 - x + 1$ will not be the same depending on the current value of $x$. If we choose a negligibility criterion of 10%, then the dominant term on $(-\infty, 10]$ and on $[10, +\infty)$ will be $x^2$, the dominant one on $[-0.1, 0.1]$ will be 1, while there will not be a single dominant term on $[-10, 10] \setminus [-0.1, 0.1]$.

We achieve this mapping by analyzing the different terms of each equation and its variables, giving the importance of the different terms in each area of the variation space. Theoretically, it is possible to extend this manipulation to logarithmic-exponential functions. However, for practical reasons, we focussed our tests on polynomials only.

First, we explore the state space by analyzing the definition space of each variable to avoid making useless computations. For example, we know for a closed chemical reaction that the total quantity of one component will not exceed the sum of the initial quantity of all the components and that its density will not exceed 1 by definition of this variable.

Using the system's equations and negligibility criteria, we automatize the model simplification with symbolic rules (particularly regarding the degree of the monomials of a polynomial equation). Polynomials can be studied using solvers to compute their zeros. Each term must be studied separately and compared to the rest of the function using the negligibility criteria and the scale chosen for this dimension of the system. These comparisons provide their ranking between known landmarks where their local orders of growth converge.

Each monomial function has an order of growth equal to its order. The resolution of the equalities between the monomial terms in a function is necessary to anticipate their dominance inversion. The symbolic rules will have to be adapted to the negligibility limit. However, we can express them as: $f(x) = o(e^{g(x)})$ for every $f$ and $g$ two polynomial functions when $x > q$, with $q$ the characteristic unitary value of the model (corresponding to the most suited magnitude for the values taken in the system), or $m(x) = o(n(x))$ if $m$ and $n$ are two monomial functions with $c(n) - c(m) = k > 0$ and $x > \sqrt{q^k}$.

In the case of multi-variable functions, we first rely on a factorization tool to separate the variable in the expression of the function. Once this step is accomplished, we treat the univariate monomials just like in a classic polynomial, treat the polynomial factors as if they were independent, and apply the simplification in the complete function. When a product of many factors cannot be simplified, we just keep it because of the many influences it merges. The result is a product of simplified polynomials that we consider our final function.

The landmarks of the partition are placed where there is equality between different components of the same equation. Around these landmarks, we place an area of qualitative equality, where even if the components are not equal, they are considered as comparable, being too close to neglect one.

Therefore, for a function $f$ composed of $n$ monomials $f_i$, we study separately the $f_i$ one after another. We inject them in a target function $f_p$ initialized at 0 such that we compute and solve the inequalities $|f_i| < |f_p|/q_f$ and $|f_i| > |f_p| * q_f$ with $q_f$ the characteristic value of $f$. In the intervals obtained by computation of these inequalities, we consider respectively that $f_i$ is negligible compared to $f_p$ (meaning that in this interval, the expression of $f_p$ will not evolve) and that $f_i$ dominates $f_p$, (we then impose that $f_p = f_i$). On the rest of the state space, we add $f_i$ to $f_p$ as they are comparable.

## C. Bounding the Error

In this section, we will present a theoretical proof that the obtained approximations give an error that we can confine or prove as negligible compared to the function itself. For zero order functions expressed as $f = g + h$ with $h$ a term or a sum of terms negligible compared to $g$ on an interval $[a, b]$, the error between $f$ and its approximation $f_p = g$ is $|h|$, that we supposed negligible. If we are interested in first order dynamics of $f$, then let us call $F = \int_a^b f(x) dx + F(a) = F_x + F(a)$, and use the same notation for $G$ and $H$. If $h$ is negligible compared to $g$ on $[a, b]$, then we can affirm that, using the intermediate value theorem, $g$ will not change its sign. As $h \ll g$, if the negligibility criteria is $k$, we have $|h| < \frac{|g|}{k}$. Therefore, $\forall t > a$, $|\int_a^t h(x) dx| < \int_a^t |h(x)| dx < \frac{\int_a^t |g(x)| dx}{k} = \frac{|\int_a^t g(x) dx|}{k}$. This means that $|H_x| < \frac{|G_x|}{k} \Rightarrow H_x \ll G_x$. This way, we can propagate the negligibility to first-order dynamics and, with a simple recurrence, to higher-order ones. Moreover, with a simple property, we can affirm that $\int_a^t |h(x)| dx < (t-a) * \max_{y \in [a,b]} (|h(y)|)$, meaning that we can bound the final error on finite rectangles. $H(a)$ must still be considered as nothing general can be expressed on it.

In the case of infinite intervals, the first part of the proof still holds. However, it is not possible anymore to limit the error as the function may not have an extremum on it.

## D. Oscillation and Random Effects Nullification

We explained a method to simplify logarithmic-exponential, and particularly polynomial functions. Our work on the order of growth and local equivalence only applies to this category of functions. However, we also propose other simplifications on other functions, taking advantage of properties of periodicity
or randomness. In the case of dynamic differential equations on long intervals, it is more complicated to assume that we can neglect a term of the equation. Even if a term is relatively negligible with respect to the main part of the equation, it can still be necessary if the simulation duration is long enough because the deviations will add up. As $\mathbb{R}$ is archimedian, adding too many terms that are negligible compared to $x^n$ may not give a result that is negligible compared to $x^n$. Consequently, reasoning dynamically on orders of magnitude on long durations has limitations. However, long-range simulations allow us to make another simplification: periodic functions can be neglected in equations if they respect some criteria regarding their amplitude and frequency. If the system is represented by an equation $f = g + h$ with $g$ non-periodic and $h$ periodic, we aim to know if it is possible to suppress or to replace $h$ in the simplified expression $f_p$ of $f$.

Let us consider the interval of study $T$. The frequency of $h$ will show if $h$ is locally periodic, i.e. if we can still observe the periodic nature of $h$ when we consider only its local behavior on $T$. When the period $p$ of $h$ is larger or comparable to $T$, then we can consider $h$ as non-periodic, and the reasoning of this paragraph does not apply. Otherwise, if $p \ll T$, the question is whether it is possible to use this periodic nature to simplify the expression. By definition, a periodic function varies in cycles around a mean value $m$, with higher and lower values that will alternate. Our logic is to consider that the lower and higher values may offset each other with enough periods in an interval. To ensure that, we must first compute $m$ on an entire period to locate its value compared to the characteristic value $q_{m,v}$ of the associated variable $v$ in the corresponding mode $m$. Let us make the hypothesis that every function we will deal with is at least piecewise continuous and, therefore, locally integrable. We can compute the medium value taken by the function on one period with $\bar{h} = \frac{\int_{t}^{t+p} h(x) dx}{p}$. If $m \ll q_v$, with $q_v$ the relative order of magnitude, then the mean effect of $h$ on the complete function can be neglected. However, we still have to look for the “variance” of $h$, represented by its amplitude. With $h_M$ and $h_m$ the extreme values of $h$ on one period, we still have to ensure that $|h_M - \bar{h}|$ and $|h_m - \bar{h}|$ are still not too important with regard to $q_v$. If $|h_M - \bar{h}| < q_v$ and $|h_m - \bar{h}| < q_v$, then our last criterion is satisfied. A high amplitude compared to $q_v$ could provoke qualitatively visible behavioral changes. If $|h_M - \bar{h}| < q_v$ and $|h_m - \bar{h}| < q_v$ but $\bar{h}$ is not negligible compared to $q_v$, we can still propose a simplification of $h$ by replacing it by its average value $\bar{h}$. If the movement around $\bar{h}$ is sufficiently limited, it is possible to consider only the mean value to identify the tendency.

When $h$ is applied to a zero-order equation, the compensation property comes easily from the periodic nature of the function and the required predicates. This means that the deviation from the mean is restricted by a value supposed at most comparable to $q_v$. In the case of higher-order differential equations, the demonstration is based on the Fourier series of periodic functions. If the function $h$ is periodic on $T$, then it can be expressed as $h(t) = \sum_{n=-\infty}^{+\infty} a_n(h)e^{-2i\pi \frac{n}{p}t}$.

Therefore, if we call $H$ the integral function of $h$, then $H = \sum_{n=-\infty}^{+\infty} -\frac{p}{2ni\pi} a_n(h)e^{-2i\pi \frac{n}{p}t}$, which means that the amplitude will regress and that the predicates that are true for zero-order functions will stay valid for their primitives if $p$ is sufficiently small.

Everything explained here for periodic functions can also be considered for stochastic terms of the functions. The analogy can be completed using the mean value of the stochastic function and by replacing the amplitude in a period with the variance of the random variable. In this case, we rely on the law of large numbers to observe the compensation after many pseudo-periods for stochastic functions. After a high number of time steps, the mean value of the random variable should converge to its theoretical mean, and the limited variance limits the potential qualitative change that could be observed. Therefore, we simplify it as earlier by replacing a stochastic function $r$, verifying these criteria by its expected value. Observing an artifact is still possible because the standard deviation does not give us a maximum amplitude but only a probabilistic value of the observed deviation to the expectancy. Nothing prevents a random variable from having an unpredictable peak value way further from the expectancy than the variance made us consider.

In brief, when faced with a periodic or stochastic term $h$ in an equation, we compute its mean value and compare it to $q_v$. If $\bar{h}$ is negligible, we consider it as 0. In any case, we must then examine the amplitude/variance of the term and its period/sampling period. If both verify the described criteria, their qualitative impact can be simplified to the mean value of the term without considering the variation around it, which will only cause local disturbances.

## IV. IMPLEMENTATION AND CASE STUDY

### A. Windy Ball and Equation Simplification

We have implemented our approach in Python, using the Sympy library to work on specific structures for polynomials. To manage periodic and stochastic functions, we manually compute the abstraction of the functions we presented earlier because of the difficulty of symbolically accessing the necessary information about the function.

To illustrate the interest of our proposition, we present two case studies on a windy ball system and a Van der Pol oscillator. The first one consists of a ball bouncing on a flat floor and exposed to the effects of an irregular wind.

We use a model where the evolution of the ball follows the ordinary differential equations:

$$ \begin{cases} \dot{x} = 10(1 + 0.5\sin(10t)) \\ \ddot{y} = -g \end{cases} $$

Starting from an initial point $(x_0, y_0) = (0, 10)$ and applying a numerical simulation from $t = 0$ to $t = 7$ with time steps of 0.007, using a second-order Runge-Kutta integration method, we get the result shown in Fig. 4. From there, abscissas and ordinates correspond to the values of $x$ and $y$.

As the term based on the $\sin$ function is periodic on $[0, 7]$ with a period $p = \frac{\pi}{5} < \frac{7}{10}$, a mean value $m = 0$ and a
<figure><img src="image_5.png" /><figcaption>Fig. 4. Windy ball reference behavior</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Windy ball : suppressing the periodic term (orange)</figcaption></figure>

maximal amplitude of 5, we can suppress this term from the
differential equation of $\dot{x}$ and compare the new result to the
first one as shown in Fig. 5. This shows that if the timing of
the impacts is not exact because of the deviation caused by
the $\sin$ term in the real system, its periodic nature makes the
deviations compensate for each other. The tendencies are well
respected, leading to a qualitative reading of the simulation
that is very close to reality.

It is possible to modify the system’s dynamics to test more
of our statements. First, we can replace the sine with a stochas-
tic term to compare the new behavior with its approximation.
In Fig. 6, we can see that with a normal term defined by a
Gaussian function $gauss(\mu, \sigma)$ with $\mu = 0.1$ and $\sigma = 1$ and
its approximation with $g(t) = \mu$, we get very close results.
This shows that it is possible to consider only the value of $\mu$
when the exposed criteria are satisfied.

We can also combine the two disturbing factors to see
whether our dominance reduction is additive. We show the
result in Fig. 7. Our approximation still efficiently captures

<figure><img src="image_7.png" /><figcaption>Fig. 6. Windy ball : suppressing the stochastic term (orange)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Windy ball : suppressing both periodic and stochastic terms</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Windy ball: hybrid system with wind shear</figcaption></figure>
the qualitative behavior of a periodic and stochastic differential equation with a simpler ODE.

Finally, we make the system more complex and transform it into a hybrid system with different equations for the wind depending on the altitude of the ball. To make this model hybrid, we suppose that the change in the wind effect is discrete and happens at a precise altitude. The equation of $\dot{y}$ does not change, but the dynamics of $\dot{x}$ becomes:

$$ \dot{x} = \begin{cases} 10(1 + 0.5\sin(10t)) & \text{if } y > 2 \\ 11 - y^2 & \text{otherwise} \end{cases} $$

The precise integration and its approximation are compared in Fig. 8. We observe a little deviation around the third bounce, but it quickly disappears and does not change the qualitative behavior. It will only change the location of some landmarks.

## B. Van der Pol and Qualitative Behavior

To highlight the interest of our proposition for qualitative behavior studies, we will apply it to the Van der Pol oscillator, which is a continuous two-dimensional physical system described by the two differential equations for the two variables $x$ and $y$, with $b$ a positive constant parameter :

$$ \begin{cases} \dot{x} = 10(y + x - \frac{x^3}{3}) \\ \dot{y} = b - x - \frac{3y}{4} \end{cases} $$

In this system, the values may oscillate depending on the value of $b$. As the oscillations happen around small values of $x$ and $y$, we choose the characteristic values $q = (1, 1)$ and choose a symmetric 2-logarithmic scale on each of them, on both positive and negative values, which create a negligibility criterion of $\frac{1}{2}$. Therefore, by applying our algorithm to these elements, we can propose a simplified version of the flow, which will then be expressed as:

$$ \begin{cases} \dot{x} = \begin{cases} 10(x + y) & \text{if } |x| < 0.5 \\ 10(y - \frac{x^3}{3}) & \text{if } |x| > 4 \\ 10(y + x - \frac{x^3}{3}) & \text{otherwise} \end{cases} \\ \dot{y} = \begin{cases} b & \text{if } |x| < b/2 \text{ and } |y| < b/2 \\ b - x & \text{if } |y| < b/2 \text{ and } |x| < 2b \\ b - \frac{3}{4y} & \text{if } |x| < b/2 \text{ and } |y| < 2b \\ -x - \frac{3y}{4} & \text{otherwise} \end{cases} \end{cases} $$

We then show the superposition of the curves obtained for a simulation of both systems with the same number of time steps and for two qualitative trajectories as shown in Fig. 9 and Fig. 10. The curves of the actual systems are in blue, and the curves of the approximated one are in orange.

We can see in both figures that despite an apparent inaccuracy in numerical values, the qualitative behaviors (convergence or limit cycles) are preserved, which is what we were aiming for. Moreover, switching from a 2-logarithmic scale to a 4-logarithmic scale almost completely removes any visible offset between the curves.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Cyclic execution of a Van der Pol system and its approximation</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Convergent execution of a Van der Pol system and its approximation</figcaption></figure>

## V. LIMITS AND NON-NULLIFIABLE EFFECTS

In this paragraph, we will expose the limitations of our contribution, which appear in borderline situations and for very sensitive systems for which it is impossible to apply the dominant term strategy we used before.

The first problem is the situation of borderline cases. When the dominated elements are in a border situation, being significantly inferior to the tendency but not below the negligibility frontier or periodic with an amplitude of the same magnitude as the characteristic quantities of the system, the approximation fails. Such terms may cause significant consequences. An example of this situation can be given by replacing the 0.5 factor of the sine function with a 2 in the windy ball case study. The mean of the complete periodic term stays unchanged, but not its amplitude. The maximum distances to the mean are moving from 5 to 20, which is no longer inferior to ten. The comparison of the simulation results in Fig. 11 illustrates this ambiguity: qualitatively, the behavior has the same nature, but the shift of the corresponding landmarks increases and even creates a phase opposition in the approximated behavior. Another example is when the system's parameters are near a behavior-switching condition. For example, in the Van der Pol oscillator presented earlier, when $b$ is close to the switching value between the convergent and divergent behavior, it is
<figure><img src="image_12.png" /><figcaption>Fig. 11. Limit of the approximation with 2 sin.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Limit of the approximation in the case of chaotic systems.</figcaption></figure>

possible to observe wrong trajectories.

The second problematic configuration comes from systems where any conditions or equations modification can completely invalidate or qualitatively change the system's behavior. A simple example derived from our previous model is a windy ball with an uneven floor. If we choose a sinusoidal floor, we can emphasize two inherent problems. The first is that it is only possible to apply the nullification to every periodic function with enough knowledge about the internal variables of this function. In our example, the cosine term takes as an argument the $x$ coordinate of the ball, which is itself a non-expressible function. It is impossible to simplify the composed function because we have no clue that $\cos x$ will indeed be periodic as we do not have an explicit expression of $x$. It is crucial to consider the composed functions, not only the top-level ones. The second problem is the sensitivity to minimal changes in the impact coordinates. Therefore, a sinusoidal term can completely change the behavior with a low amplitude. We show the results with and without approximation in Fig. 12.

With a wind of only $0.5 \sin(15t)$ added to the dynamics of $x$ in the system, we see that the behaviors of the two models of the same system diverge quickly after the first bounce. This result raises the importance of being aware of the model's sensitivity before applying a dominant simplification.

## VI. DETAILS OF THE ALGORITHM TO BUILD AN APPROXIMATION

Before concluding, we present the algorithm used to compute the qualitative approximation of the model of a system by considering each of the possible approximations presented earlier. We illustrate its application with a simple example.

This algorithm takes as entry a model in the form presented at the beginning of this article. In the following steps, the different modes of the same system will be treated separately.

The first stake is determining the unitary characteristic values $q$ for each variable and mode. If these values are given explicitly in the model, we can use them as they are. Otherwise, their computation will require more automatizing

work and rely on prior structural and symbolic knowledge of the system and an ability to deduce implicit information from it. For example, in a chemical system, it is not explicit that the masses will be preserved during the process. If we can access a variable's maximum and minimum expected values, it will be possible to deduce $q$ from them. If they have a comparable magnitude, we can easily identify the unitary value with the difference between the extrema. Otherwise, we place the unitary value at the magnitude of the minimum.

For each mode, for each equation $e$, we create a tendency function $\hat{e}$ that approximates its behavior, initialized to $\hat{e}_0 = 0$. For each known monomial term $m_i$ of $e$, we compare the magnitude order of $m_i$ with those of the terms in $\hat{e}_{i-1}$. We compare the current monomial to each term of $\hat{e}$ using the continuity of the functions by solving $\hat{e}_{i-1_j} = m_i$, $\hat{e}_{i-1_j} = m_i/k$ and $\hat{e}_{i-1_j} = km_i$, with $\hat{e}_{i-1_j}$ the $j^{th}$ term of $\hat{e}_{i-1}$ and $k$ the negligibility criteria as defined earlier. This gives us the thresholds where this term becomes negligible or dominant.

With continuity properties, we have now the position of the new monomial in comparison with each of the previous ones, and we adapt the expression of $\hat{e}$ with new segments to neglect $m_i$ when there exist $j$ with $\hat{e}_{i-1_j}$ such that $\hat{e}_{i-1_j}/m_i > k$, and to neglect every $\hat{e}_{i-1_j}$ when $m_i/\hat{e}_{i-1_j} > k$.

Theoretically, we may handle logarithmic and exponential terms by replacing $k$ with $\log k$ in exponential terms and by $e^k$ in logarithmic terms. However, this has not been tested yet.

Next, the computation continues on periodic terms of period $p$. Their mean value is computed with the formula $\bar{f} = \frac{1}{p} \int_{x}^{x+p} f(t) dt$. We then compute the amplitude with $max(\bar{f} - min_{t \in [x, x+p]}(f(t)), max_{t \in [x, x+p]}(f(t)) - \bar{f})$. Using the criteria defined in subsection III-D, we compare these values with the system unitary ones. This gives us information on the possibility of simplifying the expression of the function. Finally, this step is repeated on the stochastic terms of the equation, yielding an approximated equation with bounded expressions for specific state space areas.

For example, let us suppose that we have a system con-
sisting of 2 variables $x$ and $y$ with the equations $\dot{x} = t^3 - t^2 + 1 + y^2 - 2\sin(10t)$ and $\dot{y} = y - 2xy$. The expression of $\dot{y}$ cannot be simplified because we have no prior information about the values taken by the variables $x$ and $y$ nor any expression to compute them at any instant $t$. In the expression of $\dot{x}$, the $y^2$ term is also ignored because we do not have an explicit expression for it. Secondly, we compute the three monomial terms $t^3, t^2$ and $1$ one after another. We place ourselves in a system where the unitary value would be 1, and the negligibility criteria of 0.1. The tendency equation is initialized to 0. The first term added to $\hat{e}$ is 1. We then add the monomial $-t^2$ to $\hat{e}$. $t^2$ is negligible compared to $1$ if $|t| < \frac{1}{\sqrt{10}}$, and $1$ is negligible compared to $t^2$ if $|t| > \sqrt{10}$. After that, we can add the $t^3$ term. This one is negligible before $1$ when $|t| < \frac{1}{\sqrt{10}}$. $1$ is negligible compared to $t^3$ when $|t| > \sqrt{10}$ and $t^2$ is negligible compared to $t^3$ when $t > 10$. We then deduce a partial expression of $\hat{e}$:

$$ \hat{e} = \begin{cases} 1 + y^2 & \text{if } |t| < \frac{1}{\sqrt{10}} \\ 1 + y^2 - t^2 & \text{if } \frac{1}{\sqrt{10}} < |t| \le \frac{1}{\sqrt{10}} \\ t^3 + 1 + y^2 - t^2 & \text{if } \frac{1}{\sqrt{10}} \le |t| < \sqrt{10} \\ t^3 + y^2 - t^2 & \text{if } \sqrt{10} \le |t| \le 10 \\ t^3 + y^2 & \text{if } |x| > 10 \end{cases} $$

Finally, using the amplitude and frequency of the sine function, we can neglect $\sin(10t)$ when the absolute values taken by the other terms are superior to 2. Therefore, we add this term in the equation when $|t| < 2$.

## VII. CONCLUSION

In this paper, we have presented a methodology to reason about tendencies and dominant functions in the modeling of the behavior of complex systems:

* we analyse the ODEs to determine what are the dominant terms in different areas of the state space,

* we analyze the behavior of periodic terms in the observation window in order to determine if they can be ignored,

* we also recognize overall dominant behaviors, but that are negligible in the context of the observation window.

This study allows us to build a qualitative model of the system in the form of a map that tells us the dominant terms of the system's behavior in different areas of its state space. This qualitative map allows us to adapt the computational effort needed to simulate the system according to the position of the state of the system in the areas of the qualitative state space while preserving the qualitative behavior obtained through the simulation. This is similar to using a topographic map when hiking: if one knows they are close to a cliff, they will stop running and watch their steps, while they can progress more freely when far from any obstacle. Similarly to the hiking map, the qualitative model of the system is built only once and can be reused to optimize computational resources in simulations.

We can use the equations of the system to compute the borders of the different areas of its qualitative state space. However, we still need tools to determine beforehand the suitability of our method to a given system (depending on its

sensitivity, for instance) and the quality of the approximated qualitative behavior. This is currently the main limitation of our approach. Our long-term objective is to use this approach in the context of qualitative models design by combining it with knowledge propagation methods.

## REFERENCES

[1] M. Mitchell, Complexity: A guided tour. Oxford university press, 2009.

[2] T. A. Henzinger, "The theory of hybrid automata," in Verification of digital and hybrid systems. Springer, 2000, pp. 265–292.

[3] A. L. Brown, Qualitative knowledge, causal reasoning, and the localization of failures. MIT Artificial Intelligence Laboratory, 1974.

[4] J. De Kleer, Qualitative and Quantitative Knowledge in Classical Mechanics, ser. AI-TR-. Massachusetts Institute of Technology, Artificial Intelligence Laboratory, 1975.

[5] P. J. Hayes, "The naive physics manifesto," Expert systems in the microelectronic age, 1979.

[6] —, Readings in qualitative reasoning about physical systems. Morgan Kaufmann, 1985, ch. The second naive physics manifesto, pp. 46–63.

[7] K. D. Forbus, "Qualitative process theory," Artificial intelligence, vol. 24, no. 1-3, pp. 85–168, 1984.

[8] P. J. Mosterman and G. Biswas, "Diagnosis of continuous valued systems in transient operating regions," IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, vol. 29, no. 6, pp. 554–565, 1999.

[9] R. R. Rosenberg and D. C. Karnopp, Introduction to physical system dynamics. McGraw-Hill, Inc., 1983.

[10] B. Kuipers, "Qualitative simulation," Artificial intelligence, vol. 29, no. 3, pp. 289–338, 1986.

[11] D. Berleant and B. J. Kuipers, "Qualitative and quantitative simulation: bridging the gap," Artificial Intelligence, vol. 95, no. 2, pp. 215–255, 1997.

[12] X. Chen, E. Abraham, and S. Sankaranarayanan, "Taylor model flowpipe construction for non-linear hybrid systems," in 2012 IEEE 33rd Real-Time Systems Symposium. IEEE, 2012, pp. 183–192.

[13] J. Jerray, "Orbitador: A tool to analyze the stability of periodical dynamical systems," in ARCH@ ADHS, 2021, pp. 176–183.

[14] A. Tiwari and G. Khanna, "Series of abstractions for hybrid automata," in International Workshop on Hybrid Systems: Computation and Control. Springer, 2002, pp. 465–478.

[15] H. Zaatiti, L. Ye, P. Dague, J.-P. Gallois, and L. Travé-Massuyès, "Abstractions refinement for hybrid systems diagnosability analysis," in Diagnosability, Security and Safety of Hybrid Dynamic and Cyber-Physical Systems. Springer, 2018, pp. 279–318.

[16] S. Medimegh, "Analyse formelle de spécifications hybrides à partir de modèles sysml pour la validation fonctionnelle des systèmes embarqués," Ph.D. dissertation, Université Paris Saclay (COmUE), 2018.

[17] J.-P. Gallois and J.-Y. Pierron, "Qualitative simulation and validation of complex hybrid systems," in 8th European Congress on Embedded Real Time Software and Systems (ERTS 2016), 2016.

[18] M. L. Mavrovouniotis and G. Stephanopoulos, "Formal order-of-magnitude reasoning in process engineering," Computers & Chemical Engineering, vol. 12, no. 9-10, pp. 867–880, 1988.

[19] N. Piera and L. Travé, "About qualitative equality: Axioms and properties," in 9th International Workshop on Expert System and their Applications. Avignon, 1989.

[20] O. Raiman, "Order of magnitude reasoning," in Readings in qualitative reasoning about physical systems. Elsevier, 1990, pp. 318–322.

[21] N. Piera, M. Sanchez, and L. Travé-Massuyès, "Qualitative operators for order-of-magnitude calculers: Robustness and precision," in Thirteenth IMACS World Congress on Computation and Applied Mathematics, vol. 22, 1991.

[22] P. Dague, "Numeric reasoning with relative orders of magnitude." in AAAI, 1993, pp. 541–547.

[23] P. D. Louise Travé-Massuyès, Modèles et raisonnements qualitatifs. Hermes, 10 2003.

[24] É. Borel and A. Denjoy, Leçons sur la théorie de la croissance. Gauthier-Villars, 1910, vol. 12.

[25] G. H. Hardy, "Properties of logarithmico-exponential functions," Proceedings of the London Mathematical Society, vol. 2, no. 1, pp. 54–90, 1912.