{'fraction_non_alphanumeric': 0.05917483117998361, 'fraction_numerical': 0.016784417814903002, 'mean_word_length': 3.829406220546654, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 10, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computation of the Transient in Max-Plus Linear Systems via SMT-Solving

This paper proposes a new approach, grounded in Satisfiability Modulo Theories (SMT), to study the transient of a Max-Plus Linear (MPL) system, that is the number of steps leading to its periodic regime. Differently from state-of-the-art techniques, our approach allows the analysis of periodic behaviors for subsets of initial states, as well as the characterization of sets of initial states exhibiting the same specific periodic behavior and transient. Our experiments show that the proposed technique dramatically outperforms state-of-the-art methods based on max-plus algebra computations for systems of large dimensions.


Introduction
Max-Plus Linear (MPL) systems are a class of discrete-event systems (DES) that are based on the max-plus algebra, an algebraic system using the two operations of maximisation and addition. MPL systems are employed to model applications with features of synchronization without concurrency, and as such are widely used for applications in transportation networks [3], manufacturing [13] and biological systems [5,9]. In MPL models, the states correspond to time instances related to discrete events.
A fundamental and well-studied property of MPL systems is related to the periodic behavior of its states: from an initial vector, the trajectories of an MPL system are eventually periodic (in max-plus algebraic sense) starting from a specific event index called the transient, and with a specific period called cyclicity [3]. As explained in [13,Section 3.1], the transient is closely related to the notion of cycle-time vector, which governs the asymptotic behaviour of MPL systems.
The transient is key to solve a number of fundamental problems of MPL systems, such as reachability analysis [16] and bounded model checking [17]: it plays a crucial role as the "completeness threshold" (namely, the maximum iteration that is sufficient for the termination of the algorithm) [8] for those two problems. The computation of the transient is an interesting problem, as it is in general not correlated to the dimension of the MPL system. Thus, it is possible for the resulting transient to be relatively large for a small-dimensional MPL system. There are several known upper bounds [7,15,18,19] for the transient, which are mostly computed via the corresponding precedence graph and are, in practice, much larger than the actual values. This paper has two specific contributions. The first is to provide a novel procedure to compute the transient by means of Satisfiability Modulo Theory (SMT) solving [4]. The main idea underpinning the new method is to transform the problem instance into a formula in difference logic, and then passing the formula into an SMT solver, which outputs the transient. More precisely, in order to check the validity of the formula, we check the unsatisfiability of its negation. If the SMT solver reports "satisfied", then the original formula admits a counterexample, from which we can refine the formula. On the other hand, if SMT solver reports "unsatisfied", then from the formula we obtain the transient and the corresponding cyclicity. The second contribution of this work is to provide a procedure to synthesize the subset of the state space of an MPL system that corresponds to a specific transient/cyclicity pair. We show that one can partition the state space into sets corresponding to different transient/cyclicity pairs. The rest of the paper is structured as follows. Section 2 describes the basics of MPL systems, including the key notion of cycle-time vector. In Section 3, we provide the formal definition of transient over MPL systems and also a standard linear algebra procedure, based on matrix multiplication, to compute the transient (cf. Algorithm 1), which is later used as a benchmark. Section 4 is divided into four parts. The first part provides the background on SMT and including the underlying relevant theory. The translation of inequalities over max-plus algebra to formulae in difference logic is explained in the second part. In the third part, we provide SMT-based methods (cf. Algorithms 2 and 3) to compute the transient. The spatial synthesis problem is discussed in the last part. The comparison of the performance of the novel algorithm against the standard linear algebra procedure is presented in Section 5. The paper is concluded with Section 6. The developed code and generated data can be found in https://es.fbk.eu/people/amicheli/resources/formats20/.

Max-Plus Linear Systems
Max-plus algebra is a modification of linear algebra derived over the max-plus semiring (R max , ⊕, ⊗) where R max := R ∪ {ε := −∞} and a ⊕ b := max{a, b}, a ⊗ b := a + b, for all a, b ∈ R max . The zero and unit elements of R max are ε and 0, respectively. The max-plus algebraic operations can be extended to matrices and vectors in a natural way. For A, B ∈ R m×n max , C ∈ R n×p max and α ∈ R max , Given A ∈ R n×n max and t ∈ N, A ⊗t denotes A ⊗ . . . ⊗ A (t times). For t = 0, A ⊗0 is an n-dimensional max-plus identity matrix where all diagonal and non-diagonal elements are 0 and ε, respectively.
Given V = {v 1 , . . . , v p } as a set of vectors in R n max , we use the same notation to denote a matrix where all columns are in V i.e., The set of all max-plus linear combinations of V is called max-plus cone 1 and is denoted by cone(V ) [6]. It is formally expressed as Furthermore, we denote as v 1 , . . . , v p the basis of cone(V ). Notice that the maxplus cone is closed under the operations ⊕ and ⊗: if v, w are in cone(V ), then so is α ⊗ v ⊕ β ⊗ w for α, β ∈ R. Max-plus cones are the analogues of vector subspaces in classical linear algebra.
A dynamical system over the max-plus algebra is called a Max-Plus Linear (MPL) system and is defined as where A ∈ R n×n max is the system matrix, and vector x(k) = [x 1 (k) . . . x n (k)] encodes the state variables [3]. For example, x can be used to represent the time stamps associated to the discrete events, while k corresponds to the events counter. Hence, it is more convenient to consider R n (instead of R n max ) as the state space. Applications of MPL systems are significantly found on systems where the time variable is essential, such as transportation networks [13], scheduling or [2] manufacturing [14] problems, or biological systems [5,9].
Definition 1 (Precedence Graph [3]). The precedence graph of A ∈ R n×n max , denoted by G(A), is a weighted directed graph with nodes 1, . . . , n and an edge from j to i with weight A(i, j) for each A(i, j) = ε.
Definition 2 (Regular Matrix [13]). A matrix A ∈ R n×n max is called regular if A contains at least one finite element in each row.
Recall that a directed graph is strongly connected if, for any two different nodes i, j, there exists a path from i to j. The weight of a path p = i 1 i 2 . . . i k is equal to the sum of the edge weights in p. A circuit, namely a path that begins and ends at the same node, is called critical if it has maximum average weight, which is the weight divided by the length of the path [3].
Each irreducible matrix A ∈ R n×n max admits a unique max-plus eigenvalue λ ∈ R and a corresponding max-plus eigenspace The scalar λ is equal to the average weight of critical circuits in G(A), and E(A) can be computed from A + λ = n k=1 ((−λ) ⊗ A) ⊗k . More specifically, E(A) is the max-plus linear combination of the i th column of A + λ , for i such that A + λ (i, i) = 0 [3]. Thus, the eigenspace E(A) is a max-plus cone. A reducible matrix may have multiple eigenvalues, where the maximum one equals to the average weight of critical circuits of G(A).
which represents a simple railway network between two cities [3, Sec. 0.1], as shown in Figure 1. The dynamics w.r.t. (2) can be expressed as For 1 ≤ i, j ≤ 2, the element A(i, j) corresponds to the time taken to travel from station S j to S i , while x i (k) is the time of the k-th departure at station S i . Leaving the details aside, the matrix A has eigenvalue λ = 4 and eigenspace

Cycle-Time Vector
This section presents the definition of cycle-time vector of MPL systems. The computation of the cycle-time vector is indeed important, as it can shed light on the asymptotic behavior of MPL systems. In this section, we show its relationship with the eigenspace and eigenvalue of the underlying state matrix. Furthermore, as it will be clear in Section 3, the cycle-time vector can be used to determine whether the states of an MPL system are eventually periodic.
Definition 4 (Cycle-Time Vector [13]). Consider a regular MPL system (2), and assume that for all j ∈ {1, . . . , n} the quantity η j , defined by exists. Then the vector χ = [η 1 . . . η n ] is called the the cycle-time vector of the given sequence x(k) with respect to A.
It has been shown in [13,Theorem 3.11] that if the cycle-time vector of A exists for at least one initial vector then it exists for any initial vector. Instead of computing the limit as in Definition 4, the cycle-time vector can be generated using a procedure [11,Algorithm 31].
Theorem 1 ( [11]). Suppose we have a regular MPL system (2). For each x(0) ∈ R n there exist natural numbers p, q such that x(k + q) = (q × χ) + x(k) for all k ≥ p, where χ = [η 1 . . . η n ] is the cycle-time vector of A and the multiplication q × χ is defined in the classical algebra.

Transient in Max-Plus-Linear Systems
The transient of MPL systems is related to the sequence of the powers of matrix A, namely A ⊗k for k ≥ 0.
Proposition 2 (Transient [3,13]). For an irreducible matrix A ∈ R n×n max and its max-plus eigenvalue λ ∈ R, there exist k 0 , c ∈ N 0 , such that A ⊗(k+c) = (λ × c) ⊗ A ⊗k for all k ≥ k 0 . The smallest such k 0 and c are called the transient and the cyclicity of A, respectively.
For the rest of this paper, we denote the transient and the cyclicity of A as tr(A) and cyc(A), respectively. While cyc(A) is related to critical circuits in the precedence graph G(A) (see [3,Definition 3.94] for more details 3 ), tr(A) is unrelated to the dimension of A. Even for a small n, the transient of A ∈ R n×n max can be large. Upper bounds of the transient have been discussed in [7,15,18,19].
By Proposition 2, each irreducible MPL system enjoys a periodic behaviour with a rate λ: for each initial vector where the vectors x(1), x(2) are computed recursively by (2). A similar condition may be found on reducible MPL systems: we denote the corresponding transient and cyclicity as global, as per Proposition 2. The local transient and cyclicity for a specific initial vector x ∈ R n and for a set X ⊆ R n is defined as follows.
Definition 5. Given A ∈ R n×n max with maximum eigenvalue λ and an initial vector x ∈ R n , the local transient and cyclicity of where lcm stands for the "least common multiple".
By definition, we have tr(A, R n ) = tr(A). For a max-plus cone X = cone(V ), we show that the local cyclicity and transient can be computed from the corresponding bases, provided that tr(A, v i ) exists for all 1 ≤ i ≤ p.

Proposition 3. Given a max-plus cone
Proof. Suppose x ∈ X. It follows that there exist scalars α 1 , . . . , α p such that Definition 6. Suppose we have a regular matrix A ∈ R n×n max . The underlying MPL system (2) is classified into three categories as follows: i. never periodic: tr(A, x(0)) does not exist for all x(0) ∈ R n , ii. boundedly periodic: tr(A, x(0)) exists for all x(0) ∈ R n and tr(A) exists, iii. unboundedly periodic: tr(A, x(0)) exists for all x(0) ∈ R n but tr(A) does not.
We call (2) periodic if it is either unboundedly periodic or boundedly periodic.
We show that the periodic behavior of an MPL system is indeed related to the eigenspace and cycle-time vector of its corresponding state matrix.
Theorem 2. Suppose we have a regular matrix A ∈ R n×n max with a maximum eigenvalue λ and cycle-time vector χ. The following statements are equivalent.
Proposition 4. Suppose we have a regular matrix A ∈ R n×n max with maximum eigenvalue λ and non-empty eigenspace E(A). If there exist i ∈ {1, . . . , n} and natural numbers k 0 , c such that A ⊗k+c (·, i) = µc ⊗ A ⊗k (·, i) for all k ≥ k 0 with µ < λ then (2) is unboundedly periodic.
Proof. Because E(A) = ∅, by Theorem 2, the corresponding MPL system (2) is We now will provide the procedure to compute the transient of MPL systems. As per Proposition 2, the common method to obtain the (global) transient of A ∈ R n×n max is by computing the power of the matrix A ⊗0 , A ⊗1 , . . . until we find k 0 ≥ 0 such that A ⊗(k0+c) = λ ⊗c ⊗ A ⊗k0 where λ, c is respectively the maxplus eigenvalue and cyclicity of A. Similarly, to find the transient of A w.r.t. a max-plus cone X = cone(V ) one needs to compute A ⊗0 ⊗ V, A ⊗1 ⊗ V, . . ..
Algorithm 1 illustrates the procedure to compute transient (and cyclicity) for a max-plus cone cone(V ) w.r.t. A ∈ R n×n max . While originally designed for irreducible matrices, it also can be applied to find the transient of reducible matrices (if any). For this reason, we assign a maximum bound as termination condition. It is important to note that Algorithm 1 can also be used to compute the local transient and cyclicity for a vector: that is, when V has only one column. The algorithm starts by computing the cycle-time vector χ of the state matrix. If the entries of χ are not all the same then the transient for cone(V ) does not exist. In line 11, we perform equality checking w.r.t. a scalar between A ⊗it−m ⊗ V and A ⊗it ⊗ V . By Theorem 2 and Proposition 4, one can classify an MPL system (2) into a category in Definition 6. As a result, determining the existence of global transient is a decidable problem. For boundedly periodic MPL systems, computing the global transient is also a decidable problem. This is because they ensure the existence of a finite transient, meaning that Algorithm 1 eventually terminates. However, Algorithm 1 is sound but does not necessarily terminate (in general) for unboundedly periodic MPL systems.

Algorithm 1 Computation of cyclicity and transient of
M ← EmptyVector() empty vector used to store 3: it ← 0 number of iterations 5: χ ← CycleTimeVector(A) computing cycle-time vector 6: if elements of χ are all equal then 7: it ← it + 1 10: return it − m, m 13: if (it > N ) then 14: print "terminated after reaching maximum bound" 15: else 16: print "the transient does not exist" Remark 1. The procedure in Algorithm 1 only employs matrix operations in max-plus algebra. It can be improved by computing the cyclicity of the matrix from the corresponding precedence graph. If the resulting cyclicity is c then the range in line 10 of Algorithm 1 can be taken between 1 and c.
Example 2. Suppose we have a regular and reducible MPL system and where a ≥ 8. The corresponding eigenvalue for B is λ = 9 if 8 ≤ a ≤ 9; λ = a otherwise. Taking the power of the matrix, we have where b = max{21, 2a + 3}. One can check that, for a > 9, the matrix does not admit an eigenvector over R 3 (but it still has eigenvector over R 3 max ). As a result, B is never periodic.
On the other hand, for 8 ≤ a ≤ 9, the corresponding E(B) is not empty. Thus, B is periodic. Furthermore, for k ≥ 2, we have

From Max-Plus Algebra to Difference Logic
Before providing the main contribution, we show that the inequalities in maxplus algebra can be expressed as a formula in difference logic. For the rest of this paper, ∼ is either ≥ or >. We write ¬(a ∼ b) if it is not the case that a ∼ b.
Proposition 6 (Reduced Formula). Given real valued variables x 1 , . . . , x p and a 1 , . . . , a p , b 1 , . . . , b p ∈ R max , the inequality is equivalent to where Proof. Suppose we set initially S 1 = S 2 = {1, . . . , p}. By applying (4), (10) can be expressed as Indeed we can ignore any scalar k if a k = ε. Furthermore, for each l such that ¬(a l ∼ b l ), by (9), we have This shows that l can be removed from S 1 . Similarly, by (5), we have Again we can ignore any scalar k when b k = ε. By (6), for each l ∈ S 2 such that a l ∼ b l we have Hence l can be removed from S 2 .
Proposition 6 ensures that any inequality expression in max-plus algebra can be reduced to a simpler one in which no a variable appears on both sides i.e., S 1 ∩S 2 = ∅. However, S 1 and S 2 cannot be both empty if there exists at least one finite scalar in both sides of (10). We call (11) as a non-trivial reduced formula if both S 1 = ∅ and S 2 = ∅. (12) If S 1 = ∅ then F * ≡ false. On the other hand, if S 2 = ∅ then F * ≡ true.

Proposition 7. Given a non-trivial reduced formula in (11), then
Proof. The proof follows directly by applying Proposition 5 on (11). Proposition 7 shows that any non-trivial formula of (11) can be expressed as a difference logic formula in disjunctive and conjunctive normal forms.

Procedure to Compute Transient of MPL Systems with SMT
We now will discuss the procedure to compute the transient of an MPL system via SMT-solving. The idea behind the SMT-based procedure is to transform the equality checking in line 11 of Algorithm 1 into a formula in difference logic. Notice that the quantity M[it] in Algorithm 1 corresponds to A ⊗it ⊗ V next, and cone(V ) can be expressed as matrix V . Thus, it can be written as where p is the number of columns of V . By denoting R = A ⊗it ⊗ V and S = (λ × m) ⊗ A ⊗it−m ⊗ V , (13) can be expressed as where r ki (resp. s ki ) is the element of R (resp. S) at row k and column i. For simplicity, we denote (14) as EqFunc(R, S). By Proposition 7, each disjunct in (14) can be expressed as a formula in difference logic.
Algorithm 2 summarizes the SMT-based version of Algorithm 1. If the corresponding eigenspace of the matrix is not empty, we set the value for transient and cyclicity respectively to k 0 = 0 and c = 1 (the smallest possible for both). Then, we generate the corresponding difference logic formula F w.r.t. (13) in line 10. To check the validity of F , we use an SMT solver to check the unsatisfiability of the negation. If it is not satisfiable then the original formula is valid, and then we obtain the transient and cyclicity from the current value of k 0 and c.
On the other hand, if it is satisfiable then there exists a counterexample falsifying formula F . We express the counterexample from a satisfying assignment of ¬F as a real-valued vector w ∈ R p (line 15). Vector v = V ⊗ w corresponds to the counterexample: its transient is greater than k 0 or its cyclicity is greater than c. The resulting transient and cyclicity of v become the updated value for (k 0 , c). This process is repeated until either the SMT solver reports "unsatisfiable" in line 12 or k 0 + c exceeds the maximum bound N . (which corresponds to the termination condition of Algorithm 1).

Algorithm 2 Computation of transient and cyclicity of
if elements of χ are all equal then 4: n ← NrRows(A) 5: p ← NrCols(V ) 6: for i ∈ {1 · · · p} do 7: x[i] ← MakeSMTRealVar( ) symbolic variables 8: k0 ← 0, c ← 1 9: while ((k0 + c) ≤ N ) do 10: model ← GetSMTModel(¬F ) 12: if model = ⊥ then formula is unsatisfiable 13: return k0, c 14: else formula is satisfiable 15: if ((k0 + c) > N ) then 21: print "terminated after reaching maximum bound" 22: else 23: print "the transient does not exist" Unlike Algorithm 1, which only works on max-plus cones, Algorithm 2 can be modified (into Algorithm 3) so that it can be applied on any set of initial conditions X ⊆ R n . Although (14) is can be translated exclusively to RDL, we can extend X as an LRA formula. In line 9 of Algorithm 3, we generate a formula F which corresponds to the equality checking between A ⊗k0 and A ⊗k0+c . If X → F is valid then we have tr(A, X) = k 0 and cyc(A, X) = c. Again, to check the validity of X → F , we check the unsatisfiability of its negation. if elements of χ are all equal then 4: n ← Row(A) number of rows of A 5: for i ∈ {1 · · · n} do 6: x[i] ← MakeSMTRealVar( ) symbolic variables 7: k0 ← 0, c ← 1 8: while (k0 + c) ≤ 1000 do 9: if model = ⊥ then formula is unsatisfiable 12: return k0, c 13: else formula is satisfiable 14: v if ((k0 + c) > N ) then 19: print "terminated after reaching maximum bound" 20: else 21: print "the transient does not exist"

A Synthesis Problem
In addition to computing the transient and cyclicity of A ∈ R n×n max w.r.t. a set of initial conditions, we show that by means of difference logic and SMT, one can synthesise sets of states corresponding to specific transient (and cylicity) defined as follows On the one hand, the computation of (16) has been discussed in [1, Section 4.2] by applying backward reachability analysis. On the other hand, to the best of the authors' knowledge, there is no approach to generate (15). The following proposition shows that both (15) and (16) can be computed symbolically by expressing them as difference logic formulae: the set (15) (resp. (16)) is not empty if and only if the corresponding formula (17) (resp. (18)) is satisfiable.
Proposition 8. Given A ∈ R n×n max with global cyclicity c and maximum eigenvalue λ, we have and where Div(q) is a set of divisors of q.
Similarly, EqFunc(A ⊗p+q , λq ⊗ A ⊗p ) ∧ ¬EqFunc(A ⊗p−1+q , λq ⊗ A ⊗p−1 ) corresponds to a periodic behavior with transient p and cyclicity at most q. The remaining conjuncts guarantee that the cyclicity cannot be smaller than q.
As both (15) and (16) can be expressed as formulae in difference logic, the problem of determining the emptiness of both sets is decidable. By definition, for never periodic MPL system, S p,q = S p = ∅ for all p, q. Furthermore, for irreducible MPL systems the emptiness of (15) and (16) ii. If p > k 0 or q is not a divisor of c then S p,q (A) = ∅, iii. If S p,q (A) is empty then so is S p+1,q (A).
Proof. Let us assume the eigenvalue for A is λ. Notice that, . This shows that x(0) is an eigenvector of A ⊗c . Thus, S 0 (A) = E(A ⊗c ). The proof for S 0,1 (A) = E(A) can obtained similarly.
i. As the global transient for A is k 0 , it is by default that S p (A) = ∅ for p > k 0 . Suppose x(0) ∈ R n such that tr(A, x(0)) = k 0 . It is straightforward that x(p) ∈ S p−k0 (A) for p ≤ k 0 . This completes the proof.
ii. If p > k 0 we have S p (A) = ∅ which implies S p,q = ∅. Suppose q is not a divisor of c and S p,q = ∅. Then there exists x(0) such that cyc(A, x(0)) = lcm(c, q) > c. This contradicts the fact that c is the global cyclicity. iii. The proof is from the fact that if x ∈ S p+1,q (A) then A ⊗ x ∈ S p,q (A).
Example 3. Let us recall the 3 × 3 MPL system in Example 2 with a = 8. From the precedence graph G(B), the global cyclicity is c = 2. Leaving details aside, for p ≥ 1, we have Thus, by Proposition 8, for p ≥ 2 we have An illustration of the above sets is depicted in Figure 2.  The solid and dashed lines represent ≥ and >, respectively.

Computational Benchmarks
We compare the performance of Algorithms 1 and 3, to compute the transient of MPL systems. The experiments for both procedures are implemented in Python.
For the SMT solver, we use Yices 2.2 [10]. The computational benchmark has been implemented on an Intel R Xeon R CPU E5-1660 v3, 16 cores, 3.0GHz each, and 16GB of RAM. For the experiments, we generate 1000 irreducible matrices of dimension n, with m finite elements in each row, where the values of the finite elements are rational numbers p q with 1 ≤ p ≤ 100 and 1 ≤ q ≤ 5. The locations of the finite elements are chosen randomly. We focus on irreducible matrices to ensure the termination of the algorithms. Algorithm 1 is initialised by setting V to be a max-plus identity matrix, while for Algorithm 3 the set of initial conditions is expressed as X ≡ true.   With regards to the running time, the matrix-multiplication algorithm is faster when the values of k 0 + c are quite small. On the other hand, the larger the value of k 0 + c, the better the performance of the SMT-based algorithm is. We argue that this is because in Algorithm 3 there may be a large increment from the current guess of transient and cyclicity to the new ones. Whereas in Algorithm 1, the next candidate of transient and cyclicity is increased by one at each iteration.
As depicted in Figure 3(d), the number of finite elements m clearly affects the running time of the algorithms. We recall that the running time of Algorithm 3 depends on the satisfaction checking of a difference logic formula in line 11. The more are the finite elements, the more likely the formula is complex, and therefore the slower is the associated running time. Interestingly, based on the outcomes of the benchmarks which are presented in the Appendix, the finite elements also affect the cross-over points, which tend to increase gradually as the number of finite elements grows larger.

Conclusions and Future Work
In this paper, we have introduced a novel, SMT-based approach to compute the transient of MPL systems: our technique encodes the problem as a sequence of satisfiability queries over formulae in difference logic, which can be solved by standard SMT solvers. We have also presented a procedure to partition the statespace of MPL systems w.r.t. a given transient and cyclicity pair. The procedure has been thoroughly tested on computational benchmarks and the results show how the SMT-based algorithm is much faster that state-of-the-art techniques to compute large values of transient and cyclicity. Furthermore, we highlight that the SMT-based method can be applied to compute the transient for any initial condition, as long as it is expressible as an LRA formula.
For future research, we are interested in exploring and developing SMT-based procedures for the general model checking of MPL systems.