Improving the modified linearly implicit quantized state system methods: step size,cycle detection,and linear approximation coefficient

Abstract

The growing intricacy of contemporary engineering systems, typically reduced to differential equations, poses a difficulty in digitally simulating them. The linearly implicit quantized state system (LIQSS) provides a different method from traditional numerical integration techniques for tackling such problems. This method is effective in large sparse stiff systems and systems with frequent discontinuities. However, this method could be further improved. First, the algorithm can step through the solution analytically or through iterations. A comparison is presented in this article. Second, the intrinsic discrete behavior of this new method can cause oscillations that lead to small unnecessary simulation steps. A prior approach was made to detect and terminate these cycles. Different detection mechanisms are examined in this article. Third, a linear approximation was used. Its enhancement is also investigated in this work. Finally, the article provides which of these modifications improved the overall performance of some systems simulations using LIQSS order one.

Keywords

Quantized state systems numerical integration stiff ordinary differential equations

1. Introduction

The modeling of complex engineering systems is often reduced to differential equations that exhibit stiffness and frequent discontinuities. Stiffness and discontinuities can be a challenge in the numerical integration of such equations by classic methods.¹

Depending on the nature of the problem, traditional solvers use explicit and implicit methods. While explicit methods are cheaper, they perform poorly in the face of stiff problems, because maintaining computational stability requires a step size reduction. This increases computational cost and can make simulations impractical or even infeasible for large-scale systems.² Traditional solvers resort to implicit techniques that also involve costly iterations and matrix inversions, with an expense that grows proportionally to the size of the system.³

Quantized state system (QSS) methods are explicit and linearly implicit algorithms that aim to handle these issues more efficiently. They replace the classic time-slicing methods with a quantization of the states, leading to an asynchronous discrete-event simulation model. QSS algorithms are asynchronous because system variables do not necessarily get updated at the same time.^4,5 Solving some stiff differential equations using the linearly implicit quantized state system (LIQSS) can produce oscillations that lead to unnecessary steps. A modification was proposed to handle this issue, and a modified solver (mLIQSS) was developed and tested by Di Pietro et al.⁶

In this article, further modifications that aim to enhance the step size, the cycle detection, and the existing linear approximation, are proposed. In addition, the simulation results of these modifications are compared. The paper is organized as follows: Section 2 presents the background concepts of QSS and mLIQSS. Section 3 explains the modifications of the step size, the detection mechanism, and the linear coefficients. In section 4, the results of the modifications in the Tyson model and the advection-diffusion reaction problem are presented. Finally, section 5 concludes the article.

2. Background

2.1. QSS

QSS methods, including LIQSS, are alternatives to traditional integration methods for solving ordinary differential equations (ODEs). These methods are specially designed to solve challenges posed by the increasing complexity of modern systems. The main idea behind QSS methods is to divide the system state space into quantized regions and represent the system state in terms of these quantized values. The QSS methods update the system status only when certain thresholds or conditions are met, instead of continuously simulating the system’s behavior. These methods offer advantages such as reduced computational cost and improved performance compared with traditional numerical integration methods. In solving the following system of $n$ ODEs in Equation (1), classic methods look for the values of all variables after a small elapsed time, while QSS methods search for the time it takes to change the value of one variable by a small quantity:

\begin{matrix} \overset{\cdot}{X} = f (X, t), \end{matrix}

(1)

where $X = [x_{1}, x_{2} . . ., x_{n}]^{T}$ is the state vector, $f {: R}^{n} \times R^{+} \to R^{n}$ is the derivative function, and $t$ is the independent variable. In classic methods, the difference between $t_{k}$ (the current time) and $t_{k + 1}$ (the next time) is called the step size. In QSS, besides the step size, the difference between $x_{i} (t_{k})$ (the current value) and $x_{i} (t_{k + 1})$ (the next value) is called the quantum $Δ_{i}$ . Depending on the type of the QSS method (explicit or implicit), a new variable $q_{i}$ is set to equal $x_{i} (t_{k})$ or $x_{i} (t_{k + 1})$ respectively. $q_{i}$ is called the quantized state of $x_{i}$ .⁷ The use of the quantized state $q_{i}$ as the future value $x_{i}$ to calculate the derivatives is what makes the method implicit. Instead of solving the system shown in Equation (1), QSS methods solve the following system in Equation (2):

\begin{matrix} \overset{\cdot}{X} = f (Q, t), \end{matrix}

(2)

where $Q = [q_{1}, q_{2} . . ., q_{n}]^{T}$ is the Quantized state vector.

QSS methods were shown to have nice stability and error-bound properties⁸ while also outperforming some classic solvers.^3,9

2.2. mLIQSS

In LIQSS, all steps involve single variable updates depending on the future direction of $x_{i}$ as shown in lines (5–9) of Algorithm 1. While this intrinsic behavior of single cheap steps is advantageous, an issue arises from it. This asynchrony in updates can lead to unintended interactions between variables, causing them to stray from the true solution and trigger subsequent steps earlier than necessary. This behavior can result in cycles and unnecessary steps in the simulation. In essence, the lack of synchronization in variable updates can introduce unintended feedback loops, disrupting the simulation dynamics and potentially leading to inaccuracies or inefficiencies in the simulation results. Addressing this issue involves cycle detection mechanisms. The following conditions in Equation (3) are used to detect cycles for order 1:

\begin{matrix} {\overset{\cdot}{x}}_{j} . d x_{j} < 0 and {\overset{\cdot}{x}}_{i} . d x_{i} < 0, \end{matrix}

(3)

where $d x_{i}$ and $d x_{j}$ are linear approximations of the future derivatives. They are given as shown in Equation (4):

\begin{matrix} d x_{i} = a_{ii} . q_{i} + a_{ij} . q_{j} + u_{ij} \\ d x_{j} = a_{jj} . q_{j} + a_{ji} . q_{i} + u_{ji} . \end{matrix}

(4)

where $a_{ij}$ is the entry ( $i$ , $j$ ) of the Jacobian of the linearized derivative function and $u_{ij}$ is an affine coefficient.

Algorithm 1

mLIQSS1.^3,6

1: find the variable

i

with the fastest dynamics at the current instant
2: compute the elapsed time

e

since the last update of variable

i

3: update its value using Taylor expansion:

x_{i} = x_{i} + {\overset{\cdot}{x}}_{i} . e

4: compute the quantum using a logarithmic quantization:

Δ_{i} = \max (Δ_{i, rel}, Δ_{i, \min})

5: if

\overset{\cdot}{x} . (a_{ii} (x_{i} + sign (\overset{\cdot}{x}) . Δ) + u_{ii}) > 0

then
6:

q_{i} = x_{i} + sign (\overset{\cdot}{x}) . Δ

7: else
8:

q_{i} = \frac{- u_{ii}}{a_{ii}}

9: end if
10: compute the next time when

q_{i} = x_{i}

11: for any variable

j

depends on

i

such that

a_{ij} . a_{ji} \neq 0

do
12: compute the prediction of derivatives of

i

and

j

using Equation (4)
13: if the conditions in Equation (3) are met: then
14: perform a Backward Euler step as shown in Equation (6)
15: for any variable

k

that depends on

j

do
16: update

x_{k}

{\overset{\cdot}{x}}_{k}

a_{kj}

, and next time of change.
17: end for
18: update

a_{jj}

using Equation (7)
19: end if
20: end for
21: for any variable

j

depends on

i

do
22:

x_{j} = x_{j} + {\overset{\cdot}{x}}_{j} . e_{j}

23:

{\overset{\cdot}{x}}_{j} = f_{j} (q, t)

24: compute the next time when

q_{j} = x_{j} or | q_{j} - x_{j} | = 2 Δ_{j}

25: update

a_{ji}

using Equation (7)
26: end for
27: update

a_{ii}

using Equation (7)

A modification to the original LIQSS algorithm was implemented in Algorithm 1.

When a cycle is detected, a simultaneous update must occur. In order to simultaneously update $q_{i}$ and $q_{j}$ as the future values of $x_{i}$ and $x_{j}$ , the following equation is used:

\begin{matrix} q_{i} = x_{i} + hd x_{i} \\ q_{j} = x_{j} + hd x_{j} . \end{matrix}

(5)

Using Equations (5) and (4) is equivalent to performing a Backward Euler step. This is formulated as follows:

\begin{matrix} Q_{ij} = {(I - h . A_{ij})}^{- 1} . (X_{ij} + h . u_{ij}), \end{matrix}

(6)

where $Q_{ij} = (\begin{matrix} q_{i} \\ q_{j} \end{matrix})$ ; $A_{ij} = (\begin{matrix} a_{ii} & a_{ij} \\ a_{ji} & a_{jj} \end{matrix})$ ; $X_{ij} = (\begin{matrix} x_{i} \\ x_{j} \end{matrix})$ ; $u_{ij} = (\begin{matrix} u_{ij} \\ u_{ji} \end{matrix})$ .

The linear approximation coefficients are updated as follows:

\begin{matrix} a_{ii} = \frac{{\overset{\cdot}{x}}_{i} - x_{i}^{-}}{q_{i} - q_{i}^{-}}; a_{ij} = \frac{{\overset{\cdot}{x}}_{i} - x_{i}^{-}}{q_{j} - q_{j}^{-}} \\ a_{jj} = \frac{{\overset{\cdot}{x}}_{j} - x_{j}^{-}}{q_{j} - q_{j}^{-}}; a_{ji} = \frac{{\overset{\cdot}{x}}_{j} - x_{j}^{-}}{q_{i} - q_{i}^{-}}, \end{matrix}

(7)

where ${\overset{\cdot}{x}}_{i}^{-}$ , ${\overset{\cdot}{x}}_{j}^{-}$ , $q_{j}^{-}$ and $q_{i}^{-}$ are old values of ${\overset{\cdot}{x}}_{i}$ , ${\overset{\cdot}{x}}_{j}$ , $q_{j}$ and $q_{i}$ respectively.⁶

The value of $h$ is computed analytically or via iterations as the maximum value that satisfies Equation (8):

\begin{matrix} | q_{i} - x_{i} | \leq Δ_{i}; | q_{j} - x_{j} | \leq Δ_{j} . \end{matrix}

(8)

where $Δ_{i}$ and $Δ_{j}$ are the quantums of the variables $x_{i}$ and $x_{j}$ respectively.

2.3. Iteration and analytic approaches

Di Pietro et al.⁶ stated that the calculation of $q_{i}$ in a single update can be done using a Backward Euler step on a linear model approximation, and that when the value $h$ is finite then the update is similar to lines (5–9) of Algorithm 1. In addition, performing a backward Euler step is faster through iterations during a simultaneous update.

2.3.1. Single update

Besides the analytic approach shown in Algorithm 1, Equation (9) can be used to solve for $q_{i}$ and $h$ in a single update, such that $| q_{i} - x_{i} | \leq Δ_{i}$ and $h > 0$ :

\begin{matrix} q_{i} = x_{i} + h . d x_{i} P, \end{matrix}

(9)

where $d x_{i} P = (a_{ii} . q_{i} + u_{ii})$ is a linear approximation of the future derivative of $x_{i}$ . Taking $q_{i}$ to the left-hand side of Equation (9) leads to Equation (10):

\begin{matrix} q_{i} = \frac{x_{i} + h . u_{ii}}{1 - h . a_{ii}}, \end{matrix}

(10)

The value $h$ is the possible future step size for the variable $x_{i}$ if $d x_{i} P$ is a good prediction of the future derivative. Therefore, a search for a large $h$ value that satisfies $| q_{i} - x_{i} | \leq Δ_{i}$ could maximize the step size. On the other hand, if $d x_{i} P$ is not a good prediction of the future derivative, then the $h$ value does not depict the future step size. In this case, the distance $d_{i} = | q_{i} - x_{i} |$ may play a more important role in the future step size. In fact, the step size is determined in line (10) or line (24) of Algorithm 1. This is done by computing the next time of change, which is the time when the variable $x_{i}$ meets its quantized state $q_{i}$ or diverges from it by two quantums. That is, the value $h$ is updated again using the following equation:

h = \min (| \frac{q_{i} - x_{i} + β}{{\overset{\cdot}{x}}_{i}} |),

where $β \in {0, 2 Δ_{i}, - 2 Δ_{i}}$ . In the case the iteration method ends up with a small distance $d_{i}$ , it may produce a small future step size.

Algorithm 2 describes the iteration approach. The algorithm starts by guessing an $h$ value equal to the final simulation time minus the current simulation time. If that guess does not work, it uses the old ${\overset{\cdot}{x}}_{i}$ to compute $h$ . If that also does not work, it decreases $h$ by a linear approximation of the dependence of $| q_{i} - x_{i} |$ on h.

Algorithm 2

The Iteration approach for a single update of q_i.

h = Final Simulation Time - current Sim Time

;

q_{i} = \frac{x_{i} + h . u_{ii}}{1 - h . a_{ii}}

| q_{i} - x_{i} | > Δ_{i}

then

h = \frac{Δ_{i}}{| {\overset{\cdot}{x}}_{i} |}

;

q_{i} = \frac{x_{i} + h . u_{ii}}{1 - h . a_{ii}}

end if
if

| q_{i} - x_{i} | > Δ_{i}

then

h = h . \frac{Δ_{i}}{| q_{i} - x_{i} |}

;

q_{i} = \frac{x_{i} + h . u_{ii}}{1 - h . a_{ii}}

end if

A different explanation of the analytic approach uses Equation (10). Adding and subtracting $x_{i}$ in the right-hand side leads to Equation (11):

q_{i} = x_{i} + \frac{h . (a_{ii} . x_{i} + u_{ii})}{1 - h . a_{ii}} = x_{i} + d_{i} (h) .

(11)

For the analytic method, in order to maximize $d_{i}$ ( $maximum = Δ_{i}$ ), the mathematical function $d_{i} (h)$ shown in Figure 1 is used. Based on the positions of the pole $\frac{1}{a_{ii}}$ and the asymptote $α = \frac{- (a_{ii} x_{i} + u_{ii})}{a_{ii}}$ of this function $d_{i} (h)$ , $d_{i}$ can be determined. If $a_{ii} < 0$ and $| α | < Δ_{i}$ , then $d_{i} = | α |$ (Figure 1(a)). In other cases, $d_{i}$ can be found by the intersection between $d_{i} (h)$ and the horizontal lines $Δ_{i}$ and $- Δ_{i}$ . That is, $d_{i} = | \mp Δ_{i} |$ (Figure 1(b)).

Figure 1.

Distance $d_{i} = | q_{i} - x_{i} |$ based on the function $d_{i} (h)$ : (a) $a_{ii} < 0$ . (b) $a_{ii} > 0$ .

The analytic approach is described in Algorithm 3. At most, the algorithm checks three inequalities. This analytic method is similar to the one presented in lines (5-9) of Algorithm 1. However, it uses a different explanation that can be easily extended to higher orders, and it does not rely on the conditions of checking the signs of the derivatives.

Algorithm 3

The Analytic approach for a single update.

Compute the asymptote as

α = \frac{- (a_{ii} . x_{i} + u_{ii})}{a_{ii}}

a_{ii} < 0

then
if

α > Δ_{i}

then

q_{i} = x_{i} + Δ_{i}

else if

α < - Δ_{i}

then

q_{i} = x_{i} - Δ_{i}

else

q_{i} = \frac{- u_{ii}}{a_{ii}}

end if
else
if

α > Δ_{i}

then

q_{i} = x_{i} - Δ_{i}

else if

α < - Δ_{i}

then

q_{i} = x_{i} + Δ_{i}

else
if

α < 0

then

q_{i} = x_{i} - Δ_{i}

else

q_{i} = x_{i} + Δ_{i}

end if
end if
end if

Using the function $d_{i} (h)$ can also show the drawbacks of the iteration approach. Figure 2 shows that the iteration method can be satisfied with a value $h_{0}$ , but this value can set $q_{i}$ very close to $x_{i}$ if the asymptote is very small. In contrast, there is a value $h_{1}$ that can set $q_{i}$ at the maximum distance $Δ_{i}$ from $x_{i}$ using the analytic method.

Figure 2.

Drawbacks of the iteration approach in a single update.

2.3.2. Simultaneous update

Similar to the single update, analytic and iteration methods are introduced for the simultaneous update. To investigate the difference between them, Equation (6) is expanded. The first row leads to the following:

q_{i} = \frac{x_{i} + h . (x_{j} . a_{ij} - x_{i} . a_{jj} + u_{ii}) + h^{2} . (a_{ij} . u_{ji} - a_{jj} . u_{ij})}{1 - h . (a_{ii} + a_{jj}) + h^{2} . (a_{ii} . a_{jj} - a_{ij} . a_{ji})} .

Adding and subtracting $x_{i}$ in the right-hand side leads to Equation (12):

\begin{matrix} q_{i} = x_{i} + d s_{i} (h) \end{matrix} .

(12)

where

\begin{matrix} d s_{i} (h) = \\ \frac{h . (a_{ii} . x_{i} + a_{ij} . x_{j} + u_{ij}) + h^{2} . (a_{ij} . (a_{ji} . x_{i} + u_{ji}) - a_{jj} . (a_{ii} . x_{i} + u_{ij}))}{1 - h . (a_{ii} + a_{jj}) + h^{2} . (a_{ii} . a_{jj} - a_{ij} . a_{ji})}, \end{matrix}

is the distance between $q_{i}$ and $x_{i}$ in the simultaneous update.

Similarly, the second row leads to

\begin{matrix} q_{j} = x_{j} + d s_{j} (h), \end{matrix}

(13)

where

\begin{matrix} d s_{j} (h) = \\ \frac{h . (a_{jj} . x_{j} + a_{ji} . x_{i} + u_{ji}) + h^{2} . (a_{ji} . (a_{ij} . x_{j} + u_{ij}) - a_{ii} . (a_{jj} . x_{j} + u_{ji}))}{1 - h . (a_{ii} + a_{jj}) + h^{2} . (a_{ii} . a_{jj} - a_{ij} . a_{ji})}, \end{matrix}

is the distance between $q_{j}$ and $x_{j}$ in the simultaneous update.

Algorithm 4 describes the iteration approach in the simultaneous update. The algorithm starts by guessing an $h$ value equal to the final simulation time minus the current simulation time. If that guess does not work, it uses the old ${\overset{\cdot}{x}}_{i}$ and the old ${\overset{\cdot}{x}}_{j}$ to compute $h_{i}$ and $h_{j}$ , and it uses $h$ as the minimum of the two values. If that also does not work, it decreases $h$ . If the constraint is not satisfied, the simultaneous update is canceled and a LIQSS step is performed.

Algorithm 4

The Iteration approach in the simultaneous update order 1.

h = Final Simulation Time - current Sim Time

;

q_{i} = x_{i} + d s_{i} (h) q_{j} = x_{j} + d s_{j} (h)

if (

| q_{i} - x_{i} | > Δ_{i}

| q_{j} - x_{j} | > Δ_{j}

) then

h_{i} = \frac{Δ_{i}}{| {\overset{\cdot}{x}}_{i} |}

;

h_{j} = \frac{Δ_{j}}{| {\overset{\cdot}{x}}_{j} |}

;

h = \min (h_{i}, h_{j})

q_{i} = x_{i} + d s_{i} (h)

;

q_{j} = x_{j} + d s_{j} (h)

end if
if (

| q_{i} - x_{i} | > Δ_{i}

| q_{j} - x_{j} | > Δ_{j}

) then

h_{i} = h_{i} . \frac{Δ_{i}}{| q_{i} - x_{i} |}

;

h_{j} = h_{j} . \frac{Δ_{j}}{| q_{j} - x_{j} |}

;

h = \min (h_{i}, h_{j})

q_{i} = x_{i} + d s_{i} (h)

;

q_{j} = x_{j} + d s_{j} (h)

end if
if (

| q_{i} - x_{i} | > Δ_{i}

| q_{j} - x_{j} | > Δ_{j}

) then
cancel the simultaneous update and perform a LIQSS step
end if

For the analytic method, in order to maximize $d_{i}$ ( $maximum = Δ_{i}$ ), the position of $\mp Δ_{i}$ in respect to the mathematical function $d s_{i} (h)$ is used as shown in Figure 3. The goal is to find the intervals of acceptable values of $h \geq 0$ such that $| d s_{i} (h) | \leq Δ_{i}$ .

Figure 3.

Different positions of the quantum $Δ_{i}$ with respect to the function $d s_{i} (h)$ in a simultaneous update: (a) Zero roots, (b) One root, (c) Two roots, (d) Three roots and (e) Four roots.

The equations $d s_{i} (h) = \mp Δ_{i}$ have up to 4 roots, since the numerator and the denominator of $d s_{i} (h)$ are quadratic in $h$ and there are two right-hand sides $\mp Δ_{i}$ . Specifically, the roots are found as follows:

Let us use $b = a_{ii} x_{i} + a_{ij} x_{j} + u_{ij}$ ; $c = a_{ij} (a_{ji} x_{i} + u_{ji}) - a_{jj} (a_{ii} x_{i} + u_{ij})$ ; $α = a_{ii} + a_{jj}$ ; and $β = a_{ii} a_{jj} - a_{ij} a_{ji}$ . Then, we have

d s_{i} (h) = \frac{h . b + h^{2} . c}{1 - h . α + h^{2} . β} .

Let us start by $d s_{i} (h) = Δ_{i}$ which leads to

\begin{matrix} h . b + h^{2} . c - Δ + h . α . Δ - h^{2} . β . Δ = 0 \Leftrightarrow - Δ \\ + h . (b + α . Δ) + h^{2} . (c - β . Δ) = 0 . \end{matrix}

This last equation can have 2 real roots. Similarly, $d s_{i} (h) = - Δ_{i}$ leads to $Δ + h . (b - α . Δ) + h^{2} . (c + β . Δ) = 0$ which can have two real roots.

Based on the number of roots, acceptable intervals ( $AccInti$ ) of $h$ in which $d s_{i} (h) \in [- Δ_{i}, Δ_{i}]$ can be determined. They are indicated by the horizontal bold lines. No roots implies that $\forall h d s_{i} (h) \in [- Δ_{i}, Δ_{i}]$ , and then $AccInti = [0, \infty]$ . One root $h_{1}$ implies that $d s_{i} (h)$ has two parts, of which one is acceptable: $AccInti = [0, h_{1}]$ . Two roots $h_{1}$ and $h_{2}$ imply that $d s_{i} (h)$ has three parts; the first and the third parts are acceptable: $AccInti = [0, h_{1}] \cup [h_{2}, \infty]$ . Three roots $h_{1}$ , $h_{2}$ and $h_{3}$ imply that $AccInti = [0, h_{1}] \cup [h_{2}, h_{3}]$ . Finally, four roots $h_{1}$ , $h_{2}$ , $h_{3}$ and $h_{4}$ imply that $AccInti = [0, h_{1}] \cup [h_{2}, h_{3}] \cup [h_{4}, \infty]$ .

Similarly, acceptable intervals for variable $j$ are found ( $AccIntj$ ). After finding the accepted intervals for both values of $h_{i}$ and $h_{j}$ , the acceptable values $h$ are determined by finding the overlap between $AccInti$ and $AccIntj$ . The algorithm picks the largest segments $[L_{i}, H_{i}]$ and $[L_{j}, H_{j}]$ of $AccInti$ and $AccIntj$ respectively, where $L_{i}$ is the lower bound and $H_{i}$ is the higher bound of the segments. Next, the algorithm checks if the limit of one segment is in the other segment. In the case there is an overlap, the algorithm picks the optimal limit, and if it is neither null nor infinite, it computes $q_{i}$ and $q_{j}$ using the functions $d s_{i} (h)$ and $d s_{j} (h)$ . Otherwise, $q_{i}$ and $q_{j}$ are picked as the asymptotes of the functions $d s_{i} (h)$ and $d s_{j} (h)$ . In the case there is no overlap, the algorithm checks the next segments in $AccInti$ and $AccIntj$ . The overlap in one of the segments is guaranteed since $d s_{i} (0) = 0$ and $d s_{j} (0) = 0$ . This approach is presented in Algorithm 5.

Algorithm 5

The analytic approach in the simultaneous update order 1.

construct the set of acceptable intervals

AccInti

construct the set of acceptable intervals

AccIntj

while

AccInti and AccInti

have segments do ▷ at most 3 iterations
Pick the largest segments

[L_{i}, H_{i}]

and

[L_{j}, H_{j}]

AccInti

and

AccIntj

H_{i} \in [L_{j}, H_{j}] ∥ H_{j} \in [L_{i}, H_{i}]

then ▷ overlap exists

h = \min (H_{i}, H_{j})

h = \infty & (Lj \neq 0 | | Li \neq 0)

then

h = \max (Lj, Li)

end if
if

h = \infty

; then ▷ case: both zero roots

q_{i} = x_{i} + asymptote of d s_{i} (h) and q_{j} = x_{j} + asymptote of d s_{j} (h)

else

q_{i} = x_{i} + d s_{i} (h) and q_{j} = x_{j} + d s_{j} (h)

end if
else
Remove the larger segment
end if
end while

A search for the largest overlap between $AccInti$ and $AccIntj$ takes at most three iterations because the maximum number of segments in $AccInti$ and $AccIntj$ is three. In addition, the analytic method guarantees that when there are roots, at least one distance $d_{i} = | q_{i} - x_{i} |$ or $d_{j} = | q_{j} - x_{j} |$ is equal to the corresponding quantum.

3. Modifications

This section presents two possible modifications to the current mLIQSS algorithm. First, different conditions to predict cycles are investigated. Second, a different approach to updating the linear coefficients is introduced.

3.1. Cycle detection mechanism

The fundamental concept of the QSS method revolves around the cheap update of a single variable. However, such an update can trigger sudden shifts in other variables, resulting in shorter steps. The reduction of the step count caused by a simultaneous update within a trajectory is a sign of the presence of a cycle. Yet, practically this criterion cannot be directly exploited as we lack the forthcoming step count. Hence, each step necessitates a judicious decision. Currently, order 1 cycle detection relies on the change of both derivatives’ signs as shown in Equation (3). These criteria may miss cycles in some cases.

First, consider the case shown in Figure 4. This situation contains cycles, but at every step, only one variable changes direction. At $time = t 1$ , $x_{i} = q_{i}^{-}$ and $q_{i}$ is set up without a change of direction. At $time = t 2$ , the increase of $q_{i}$ causes $x_{j}$ to shift down, and $d x_{j} (2)$ keeps the same sign as ${\overset{\cdot}{x}}_{j} (2)$ . At $time = t 3$ , the decrease of $q_{j}$ causes $x_{i}$ to shift down, and $d x_{i} (3)$ keeps the same sign as ${\overset{\cdot}{x}}_{i} (3)$ . At $time = t 4$ , the decrease of $q_{i}$ causes $x_{j}$ to shift up, but $d x_{j} (4)$ keeps the same sign as ${\overset{\cdot}{x}}_{j} (4)$ . A significant change with a direction change can be considered for cycle detection as shown in Equation (14):

\begin{matrix} | {\overset{\cdot}{x}}_{j} - d x_{j} | > \frac{| {\overset{\cdot}{x}}_{j} + d x_{j} |}{2} and | {\overset{\cdot}{x}}_{i} - d x_{i} | > \frac{| {\overset{\cdot}{x}}_{i} + d x_{i} |}{2} . \end{matrix}

(14)

Figure 4.

Using the change of the derivatives’ signs in cycle detection.

Second, consider the case shown in Figure 5. At $time = t 1$ , $q_{i}$ is set up from $x_{i}$ using the predicted derivative $d x_{i} P = a_{ii} q_{i} + u_{ii}$ in the single update. This update assumes that the influence of the other variable $j$ is negligible. $d x_{j}$ and $q_{j}$ shift down because of this $q_{i}$ . The estimate of the future derivative $d x_{i}$ is not significantly different from ${\overset{\cdot}{x}}_{i}$ . Therefore, the conditions for cycle detection are not met. Then at $time = t 2$ , $q_{j 2}$ is set up using $d x_{j} P = a_{jj} q_{j} + u_{jj}$ , and ${\overset{\cdot}{x}}_{j}$ shifts up. $x_{iNEW}$ shifts up because of this new value $q_{j 2}$ . In this situation, $q_{j}$ influences ${\overset{\cdot}{x}}_{i}$ , but the comparison of $d x_{i}$ and ${\overset{\cdot}{x}}_{i}$ did not capture this influence. It can be more efficient to compare $d x_{i} P$ and $d x_{i}$ since $d x_{i} P$ is the first predicted derivative that assumes the influence of the variable $j$ is negligible. Then, $d x_{i}$ verifies the validity of this assumption.¹⁰

Figure 5.

Comparing $d x_{i} P$ and $d x_{i}$ in cycle detection.

Therefore, Equation (15) can also be used for cycle detection:

\begin{matrix} {\overset{\cdot}{x}}_{j} . d x_{j} < = 0 and d x_{i} P . d x_{i} < = 0 . \end{matrix}

(15)

Finally, both former cases can be considered as shown in Equation (16):

\begin{matrix} | {\overset{\cdot}{x}}_{j} - d x_{j} | > \frac{| {\overset{\cdot}{x}}_{j} + d x_{j} |}{2} and | d x_{i} P - d x_{i} | > \frac{| d x_{i} P + d x_{i} |}{2} . \end{matrix}

(16)

3.2. More accurate linear coefficients

For accurate identification of cycles and optimal step sizes, the prediction of future derivatives through a linear approximation must be highly precise. Erroneous values could force the step size to be small, miss cycles, or conduct unnecessary expensive simultaneous updates, which could hurt the overall performance. Therefore, the coefficients of the linear approximation $a_{ii}, a_{ij}, a_{jj}, a_{ji}$ , $u_{ii}, u_{ij}$ and $u_{ji}$ need to be as accurate as possible. These coefficients are Jacobian entries of the linearization function shown in Equation (19) of the derivative function $f_{i}$ shown in Equation (17):

{\overset{\cdot}{x}}_{i} = f_{i} (q, t),

(17)

where $q$ is the vector of all the quantized variables that $x_{i}$ depends upon, and $t$ is the independent variable. The current approach of finding the linear coefficients has two drawbacks.

First, the method used in Equation (7) is not accurate in case of a simultaneous update of $q_{i}$ and $q_{j}$ . In fact, the derivative ${\overset{\cdot}{x}}_{i}$ and its old value ${\overset{\cdot}{x}}_{i}^{-}$ presented in Equation (19) and Equation (18) respectively do not lead to the expression that finds $a_{ii}$ in Equation (7), when subtraction is conducted:

{\overset{\cdot}{x}}_{i}^{-} = a_{ii} . q_{i}^{-} + a_{ij} . q_{j}^{-} + u_{ij},

(18)

{\overset{\cdot}{x}}_{i} = a_{ii} . q_{i} + a_{ij} . q_{j} + u_{ij} .

(19)

Another variable ${\overset{\cdot}{x}}_{i_{new}}^{-}$ that accounts for the change of $q_{j}$ can be used to compute $a_{ii}$ :

{\overset{\cdot}{x}}_{i_{new}}^{-} = a_{ii} . q_{i}^{-} + a_{ij} . q_{j} + u_{ij} .

(20)

Subtracting Equation (20) from Equation (19) leads to the expression that finds $a_{ii}$ in Equation (7). While this modification enhances accuracy, it is expensive as it costs an extra function call for every coefficient $a_{ii}, a_{ij}, a_{jj}$ , and $a_{ji}$ .

Second, in large nonlinear systems where the coefficient $a_{ii}$ depends on other variables that change after its update, the coefficient becomes outdated. Consider the example shown in Equation (21) and in Figure 6:

{\overset{\cdot}{x}}_{i} = f_{i} (q_{i}, q_{j}, q_{k}, t) = (1 + q_{j} + q_{k}) . q_{i} + q_{j} + u_{ij} .

(21)

Figure 6.

$a_{ii}$ is outdated.

At the end of the step at time $t_{1}$ , the linearization coefficient $a_{ii}$ is computed using the expression that finds $a_{ii}$ in Equation (7). At time $t_{2}$ , $q_{j}$ and $q_{k}$ can undergo a simultaneous change, which affects the value of ${\overset{\cdot}{x}}_{i}$ . In the next step $t_{3}$ , the variable $x_{i}$ needs an update, and it uses $a_{ii}$ , which is no longer valid.

The article investigates another possible way to approximate the derivative function using the Jacobian entries of the original derivative function as coefficients of the linearization $a_{ii} q_{i} + a_{ij} q_{j} + u_{ij}$ , via Equation (22),

a_{ii} = \frac{\partial f_{i}}{\partial q_{i}} .

(22)

Hence, $a_{ii}$ can be computed from the Jacobian expression rather than from Equation (7). By using an expression for the coefficients, we can update them at the time we need them.

Going back to the example in Equation (21), the coefficient $a_{ii}$ can be calculated as $a_{ii} = \frac{\partial f_{i}}{\partial q_{i}} = 1 + q_{j} + q_{k}$ . Each time $a_{ii}$ is needed, a simple update using this latter expression would account for the change of $q_{j}$ and $q_{k}$ .

Computing the original Jacobian entries of the system can be better and cheaper than finding the coefficients manually, as finding these values costs one function call. This approach extracts the Jacobian entries as expressions in a function during compilation time. Hence, this modification can balance accuracy with computational efficiency in the simulation of complex dynamic systems.

3.3. Implementation

A Linearly implicit algorithm is implemented in a stand-alone solver programmed in the Julia language. This work is inspired by the previous work using the stand-alone solver in the C language.¹¹ The code can be found at https://github.com/mongibellili/QSSolver.git.

4. Results

First, this section compares the iteration and the analytic approach for updating the quantized variable. Second, it presents the simulation results of the two modifications using the Tyson and Advection-Reaction-Diffusion (ADR) problems. Third, the performance of the final version is compared with the former mLIQSS algorithm in Julia and the mLIQSS in the C language. Finally, the performance of the final version is compared with the Implicit Euler from the DifferentialEquations.jl package.¹²

The comparison includes the number of steps, relative error, and CPU time using an Intel(R) Core(TM) i5-1135G7 @ 2.40 GHz under Windows 11. The relative error is computed as: $err = \frac{\sum_{i = 1}^{n} er r_{i}}{n}$ , with

er r_{i} = \sqrt{\frac{\sum_{k = 1}^{m} {(so l_{i} [k] - re f_{i} [k])}^{2}}{\sum_{k = 1}^{m} re f_{i} {[k]}^{2}}},

(23)

where $n$ is the number of variables in the system, $so l_{i}$ is the obtained solution of the $i^{th}$ variable (a vector of length $m$ ), $re f_{i}$ is the reference solution of the $i^{th}$ variable. The reference solution is obtained from the DifferentialEquations.jl package at absolute tolerance abstol = $10^{- 12}$ and relative tolerance reltol = $10^{- 8}$ using the solver Feagin14 (a Runge-Kutta method of order14) which is available at Rackauckas and Nie.¹² The notations in Table 1 are used in the different Tables presenting the results:

Table 1.

Names of different variations.

Name	Meaning
coef1	coefficients as Jacobian entries of the linearization function
coef2	coefficients as Jacobian entries of the original derivative function
single1	Single quantized variable update using an analytic approach
single2	Single quantized variable update using a second explanation of the analytic approach
single3	Single quantized variable update using iterations
detect1	cycle detection conditions using Equation (3)
detect2	Cycle detection conditions using Equation (14)
detect3	Cycle detection conditions using Equation (15)
detect4	Cycle detection conditions using Equation (16)
simul1	Simultaneous update using iterations
simul2	Simultaneous update using an analytic approach

coef1_single1_simul1_detect1 is the former mLIQSS algorithm.

4.1. The Tyson model

The first problem is a model that is of great interest in the field of biology. The Tyson model describes the mechanism of a cell: the division, the growth, and the interactions among its components.⁶ The model contains six components ( $C 2, CP, pM, M, Y, YP$ ), which are the variables of the following system of equations in Equation (24):

\begin{matrix} \frac{C 2}{dt} = k_{6} M - k_{8} C 2 + k_{9} CP \\ \frac{CP}{dt} = k_{3} CPY + k_{8} C 2 - k_{9} CP \\ \frac{pM}{dt} = k_{3} CPY - pM {(k_{1} 0 + k_{4} \frac{M}{(C 2 + CP + pM + M)})}^{2} + k_{5} M \\ \frac{M}{dt} = pM {(k_{1} 0 + k_{4} \frac{M}{(C 2 + CP + pM + M)})}^{2} + k_{5} M - k_{6} M \\ \frac{Y}{dt} = k_{1} - k_{2} Y - k_{3} CPY \\ \frac{YP}{dt} = k_{6} M - k_{7} YP, \end{matrix}

(24)

where $k_{i}$ are parameters. The model variation that is tested in this article uses the following values: $k_{1} = 0.015; k_{2} = 0; k_{3} = 200; k_{4} = 180; k_{5} = 0; k_{6} = 1; k_{7} = 0.6; k_{8} = 10^{6}; k_{9} = 10^{3}; k_{10} = 0.018$ . More details about this model can be found in Tyson.¹³

We simulated this system at tolerance values $abstol = 10^{- 3}$ and $reltol = 10^{- 2}$ until a final time $t_{f} = 25$ .

4.1.1. Single and simultaneous updates

Table 2 compares the iteration and the analytic approach. The number of steps is slightly higher in the single update with iterations due to possible small distances between $q_{i}$ and $x_{i}$ . The results look similar for the second explanation of the analytic approach. Using the analytic simultaneous update leads to a higher step count. This is because to solve for $q_{i}$ , $q_{j}$ and $h$ analytically might push $q_{i}$ or $q_{j}$ away from $x_{i}$ or $x_{j}$ by $Δ_{i}$ or $Δ_{j}$ while the system is near equilibrium. Therefore, maximizing the distance $d_{i} = | q_{i} - x_{i} |$ in a simultaneous update can cause a slower convergence near an equilibrium.

Table 2.

Comparison of the number of steps of the iteration and the analytic approaches.

Updates	Simul1	Simul2
single1	1337	1555
single2	1337	1555
single3	1376	1382

4.1.2. Cycle detection mechanism

Table 3 shows that there is no improvement in any of the different cycle detection mechanisms while using the former method of updating the linear coefficients.

Table 3.

Comparing different cycle detection conditions.

Cycle detection condition	Number of steps	Number of simultaneous steps
detect1	1337	1
detect2	1371	229
detect3	1365	91
detect4	1488	314

4.1.3. Linear approximation coefficient

Table 4 shows that the use of the Jacobian of the original derivative function is slightly better than the former method of computing the linear coefficients because the two drawbacks of the former method do not occur frequently. In this case, a simultaneous update is performed only once. On the other hand, if the detection mechanism $detect 3$ is used, more simultaneous updates are performed, resulting in better performance when using the more accurate coefficients as shown in Table 5.

Table 4.

Comparing the different calculations of the linear coefficient.

Linear coefficient	Number of steps	Number of simultaneous steps
Coef1	1337	1
Coef2	1317	1

Table 5.

Variations of the solver with coef2.

Variation	Number of steps	Number of simultaneous steps
single1 and detect1	1317	1
single1 and detect3	1287	66

4.1.4. Comparison of final version with former version and other solvers

Table 6 shows the results of the final version (coef2_single1_simul1_detect3) of mLIQSS compared with the former version as well as other state-of-the-art solvers. The two modifications improved slightly the final mLIQSS (Julia), and with close results to the final version of mLIQSS1 (c). All mLIQSS solvers did not perform better than the Implicit Euler because the Tyson model is not a large sparse system.

Table 6.

Comparison of final version.

Solver	Number of steps	CPU time (ms)	Error ( $abstol = 10^{- 5}$ )
Former version (Julia)	1337	0.794	0.0969
mLIQSS1 (c)	1024	0.726	0.0976
Final version (Julia)	1285	0.786	0.0990
Implicit Euler (Julia)	70	0.171	0.0431

4.2. The Advection Diffusion Reaction problem

The Advection Diffusion Reaction equations describe many processes that include heat transfer, chemical reactions, and many phenomena in areas of environmental sciences. They are ODEs that result from the method of lines (MOL). The resulting system in Equation (26) is a large stiff system with possible large entries outside the main diagonal:⁶

\begin{array}{l} For i = 1 \dots N - 1 {\dot{u}}_{i} = - a \frac{u_{i} - u_{i - 1}}{Δ x} + d \frac{u_{i + 1} - 2 u_{i} + u_{i - 1}}{Δ x^{2}} + r (u_{i}^{2} - u_{i}^{3}) \\ {\dot{u}}_{N} = - a \frac{u_{N} - u_{N - 1}}{Δ x} + d \frac{2 u_{N - 1} - 2 u_{N}}{Δ x^{2}} + r (u_{N}^{2} - u_{N}^{3}) \end{array}

(25)

where N is the number of grid points and $Δ x = \frac{10}{N}$ is the grid width after the discretization of the problem with the MOL, $a$ is the advection parameter, $d$ is the diffusion parameter, and $r$ is the reaction parameter. The initial condition is given by

u_{i} (t = 0) = {\begin{matrix} 1 & if i \in [1, N / 3] \\ 0 & else . \end{matrix} .

Further explanations can be found in Bergero et al.¹⁴

We simulated this system until a final time $t_{f} = 10$ at tolerance values $abstol = 10^{- 2}$ and $reltol = 10^{- 1}$ . The advection parameter is fixed at $a = 1$ , and the reaction parameter is fixed at $r = 1000$ . The number of grid points is picked as $N = 1000$ and $d$ is set to $0.1$ .

4.2.1. Single and simultaneous updates

As seen in the Tyson model, Table 7 shows that the number of steps in the iteration approach is higher for the single updates and lower for the simultaneous update. The number of steps is not exactly the same for the two analytic approaches because they differ in the way they rely on the coefficient $a_{ii}$ , which can be inaccurate in some steps.

Table 7.

Comparison of the number of steps of the iteration and the analytic approaches.

Updates	Simul1	Simul2
single1	712072	4465853
single2	708299	4455364
single3	1302669	5171775

4.2.2. Cycle detection mechanism

Table 8 shows that there is no improvement in any of the different cycle detection mechanisms while using the former method of computing the linear coefficients.

Table 8.

Comparing different cycle detection conditions.

Cycle detection condition	Number of steps	Number of simultaneous steps
detect1	712072	8304
detect2	7056848	3527398
detect3	2696186	757380
detect4	7634858	4274875

4.2.3. Linear approximation coefficient

Table 9 shows that the benefit of using the Jacobian of the original derivative function is more obvious in this problem as the number of variables and the number of simultaneous updates are higher. The final version is better because it uses more accurate coefficients while performing more simultaneous updates. In addition, the results of the two analytic approaches are the same because the linear coefficient is more accurate, and the cycle detection mechanism $detect 3$ is effective in this final version of the solver as shown in Table 10.

Table 9.

Comparing the different calculations of the linear coefficient.

Linear coefficient	Number of steps	Number of simultaneous steps
Coef1	712072	8304
Coef2	82652	1222

Table 10.

Variations of the solver with coef2.

Variation	Number of steps	Number of simultaneous steps
single1 and detect1	82652	1222
single3 and detect1	82652	1222
single1 and detect3	43224	13252

4.2.4. Comparison of final version with former version and other solvers

Table 11 shows the results of the final version of mLIQSS compared with the former version as well as other state-of-the-art solvers. First, this problem shows clearly the advantages of the two modifications in the final version of the solver. The final version of the solver is better than the former version and the mLIQSS1 (c) solver. Two different CPU times are recorded for this latter solver. The first time is the total run time including the initialization, simulation, and save times. The second is only the simulation time. In addition, all mLIQSS solvers are better than the Implicit Euler because the Advection Diffusion Reaction problem is a large sparse system.

Table 11.

Comparison of the final version.

Solver	Number of steps	CPU time (ms)	Error
Former version (Julia)	712072	1548.42	0.1369
mLIQSS1 (c)	107073	645.272 (56.412)	0.2290
Final version (Julia)	43224	94.99	0.0268
Implicit Euler (Julia)	481	13134	0.0572

5. Conclusion

Quantization-based methods demonstrated encouraging results, particularly in dealing with large sparse stiff models and scenarios characterized by frequent discontinuities. Despite considerable research that aimed to enhance their efficiency, there remains scope for further improvement, as quantization-based approaches continue to be an evolving area of study. In this work, we performed the following tasks. First, an investigation of current approaches to update the quantized variable is conducted. Then, two modifications are introduced and tested in the Tyson model and the Advection Diffusion Reaction system. The main findings are as follows:

An iteration approach in the single update of the quantized state performed worse than the existing analytic approach because the iteration approach can lead to small distances between the variable $x_{i}$ and its quantized state $q_{i}$ .

An alternative explanation to the existing analytic approach, which can be more easily extended to higher orders, is presented and tested.

The analytic approach in the simultaneous update proved ineffective because maximizing the distance $d_{i} = | q_{i} - x_{i} |$ near an equilibrium can slow down the convergence.

Through an integration of the Jacobian entries of the original system, a more precise cycle detection and an adaptive step size are achieved.

Modifications of the cycle detection mechanism slightly improved performance.

With this investigation, QSS methods further solidify their position as powerful tools for simulating dynamic systems across various domains where traditional methods may struggle to maintain accuracy and efficiency. Regarding future work, we are considering the following:

Investigating the effect of the modifications on higher order solvers and in other stiff systems.

Conducting an error analysis considering these sources of error: the linear approximation, the quantization, and the cycle detection, in conjunction with the work presented in Mahawattege and Nutaro.⁸

Footnotes

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

ORCID iD

Elmongi Elbellili

Author biographies

Elmongi Elbellili graduated from the United States Military Academy, New York as an Electrical Engineer in 2008. Then, he obtained a Master of Science in computer science from Naval Postgraduate School California in June 2015. He is currently a PhD student in the department of computer science at KU Leuven Belgium. His research includes quantized state methods, numerical integration, differential equations, and simulation.

Philippe Blondeel obtained his PhD in Engineering Mathematics from the Department of Computer Science at KU Leuven in 2022. His PhD focused on Monte Carlo methods and their variants, such as Multilevel Monte Carlo (MLMC) and Multilevel quasi-Monte Carlo (MLQMC), with applications in uncertainty quantification. He is currently working at the Belgian Royal Military Academy as a researcher in the Mathematics department. His current research encompasses, among others topics, (quasi-)Monte Carlo methods, optimization methods and stochastic optimal control methods.

Daan Huybrechs is a research professor in Mathematical Engineering at the Department of Computer Science at KU Leuven. His research focuses on function approximation, oscillatory integrals, and the simulation of oscillatory integral equations, with applications in wave problems from acoustics and electromagnetics. He obtained his Ph.D. in 2006 and M.Sc. in 2002, both from the Department of Computer Science at KU Leuven. He is currently an associate editor of the IMA Journal of Numerical Analysis.

Ben Lauwens graduated from the Royal Military Academy, Belgium, as a telecommunications engineer in 2000. He obtained a Master of Science in Information and Communication Technologies from KU Leuven in 2005 and a PhD in engineering science from the Electrical Engineering Department at KU Leuven in 2009. He is currently a professor and head of the Department of Mathematics at the Royal Military Academy. His research includes modeling, computer simulation, optimization (including artificial intelligence), and quantum computing.

References

Hairer

Wanner

Solving ordinary differential equations II: stiff and differential-algebraic problems. Berlin: Springer, 1991.

Hairer

Nørsett

Wanner

Solving ordinary differential equations I: nonstiff problems. 2nd ed. Berlin: Springer, 1993.

Cellier

Kofman

Continuous system simulation. New York: Springer, 2006.

Migoni

Bortolotto

Kofman

, et al. Linearly implicit quantization-based integration methods for stiff ordinary differential equations. Simul Model Pract Theory 2013; 35: 118–136.

Castro

Bergonzi

Marcosig

, et al. Discrete-event simulation of continuous-time systems: evolution and state of the art of quantized state system methods. Simulation 2024; 100: 613–638.

Di Pietro

Migoni

Kofman

. Improving linearly implicit quantized state system methods. Simulation 2019; 95: 127–144.

Migoni

Kofman

Cellier

Quantization-based new integration methods for stiff odes. Simulation 2012; 88: 387–407.

Mahawattege

Nutaro

Aquantized state integrator with second order errors over monotonic segments. In: 2022 annual modeling and simulation conference (ANNSIM), San Diego, CA, 18–20 July 2022, pp. 428–436. San Diego, CA: Society for Modeling & Simulation International (SCS).

Kofman

Relative error control in quantization based integration. Lat Am Appl Res 2009; 39: 231–238.

10.

Elbellili

Lauwens

Huybrechs

Investigation of different conditions to detect cycles in linearly implicit quantized state systems. In: Proceedings of the 3rd international conference on computational modeling, simulation and optimization, Phuket, Thailand, 14–16 June 2024, pp. 469–473. New York: IEEE.

11.

Fernández

Kofman

A stand-alone quantized state system solver for continuous system simulation. Simulation 2014; 90: 782–799.

12.

Rackauckas

Nie

Differentialequations.jl—a performant and feature-rich ecosystem for solving differential equations in Julia. J Open Res Softw 2017; 5: 15–24.

13.

Tyson

Modeling the cell division cycle: cdc2 and cyclin interactions. Proc Natl Acad Sci U S A 1991; 88: 7328–7332.

14.

Bergero

Fernández

Kofman

, et al. Time discretization versus state quantization in the simulation of a 1d advection-diffusion-reaction equation. Simulation 2016; 92: 47–61.