Switched max-plus linear-dual inequalities: cycle time analysis and applications

P-time event graphs are discrete event systems suitable for modeling processes in which tasks must be executed in predefined time windows. Their dynamics can be represented by max-plus linear-dual inequalities (LDIs), i.e., systems of linear dynamical inequalities in the primal and dual operations of the max-plus algebra. We define a new class of models called switched LDIs (SLDIs), which allow to switch between different modes of operation, each corresponding to a set of LDIs, according to a sequence of modes called schedule. In this paper, we focus on the analysis of SLDIs when the considered schedule is fixed and either periodic or intermittently periodic. We show that SLDIs can model a wide range of applications including single-robot multi-product processing networks, in which every product has different processing requirements and corresponds to a specific mode of operation. Based on the analysis of SLDIs, we propose algorithms to compute: i. minimum and maximum cycle times for these processes, improving the time complexity of other existing approaches; ii. a complete trajectory of the robot including start-up and shut-down transients.


Introduction
Time-window constraints arise in different settings. Typical examples from manufacturing are the bakery and the semi-conductor industries, where some processes (e.g., yeast fermentation in the former case, metal deposition in the latter) must be executed under strict temporal constraints in order to obtain the desired quality of the final product [1,2]. In transportation, time windows can be used to specify admissible pickup and delivery times for customers [3]. From the pharmaceutical domain, high-throughput screening systems are another example where, due to chemical requirements, the time allowed to elapse between two operations is restricted by lower and/or upper bounds [4]. P-time event graphs (P-TEGs) are event graphs 1 in which tokens are forced to sojourn in places in predefined time windows. Given their ability to model time-window constraints, they have been applied to solve scheduling and control problems in various types of systems including bakeries, electroplating lines, and cluster tools [5][6][7][8]. The signal describing firing times of transitions in P-TEGs evolves -non-deterministically -according to max-plus linear-dual inequalities (LDIs), i.e., dynamical inequalities that are linear in the primal operations and in the dual operations of the max-plus algebra (see (2) in Section 2.4).
Since temporal upper bound constraints can be considered as specifications that a system driven exclusively by lower time bounds (i.e., a max-plus linear system) needs to satisfy, several authors have addressed the problem of limiting the sojourn time of tokens in places from a control point of view. This resulted in a variety of techniques; we mention for instance [9][10][11][12], in which the authors find sufficient conditions for the control of max-plus linear systems subject to constraints using, respectively, geometric control theory, max-plus spectral theory, residuation theory, and model predictive control. Other sufficient conditions were discovered in [13]. Note that, in its most general formulation, this control problem is still open, as no necessary and sufficient condition for its solution has been found. The reason for this can be traced to the fact that a more fundamental problem, namely, checking the existence of solutions of LDIs, has not yet been fully understood. On the other hand, the steady-state version of the problem was solved for a large class of instances in [14], exploiting the results of [9]. Moreover, the simple nature of the dynamics of P-TEGs has led to a number of theoretical results regarding the analysis of their cycle time [15][16][17][18], which is defined as the temporal difference between the occurrence of two repetitions of the same event in the system, assuming events occur in a periodic manner (for a formal definition, see Section 3.2).
Such a simple characterization comes, however, at the cost of limited modeling power. To illustrate this, consider the example of flow shops, which are manufacturing systems consisting of a sequence of machines of unitary capacity where all the jobs (i.e., parts to be processed) must visit all the machines in the same order (see Figure 1). P-TEGs are the ideal tool for representing flow shops where all the jobs are identical. Suppose instead that jobs of different types are present, each of which requires different processing times in machines, and that the entrance order of jobs is periodic and repeats every V ∈ N jobs. An example of periodic entrance order of period V = 2 for the flow shop of Figure 1 is ababab . . . In this case, P-TEGs can still be adopted to model the system (this is shown formally in Proposition 5), but: • the number of transitions in the resulting P-TEGs increases (linearly) with V , which leads to a rather high computational complexity for the cycle time analysis where large values of V are considered; • for a different entrance order of jobs, a new P-TEG must be built. Moreover, when the number of jobs to be processed is infinite 2 and the entrance order is not periodic, no P-TEG can model the manufacturing system.
With the aim of overcoming these limitations, in this paper we introduce a new class of dynamical systems called switched max-plus linear-dual inequalities (SLDIs). SLDIs extend the modeling power of P-TEGs by allowing to switch among different modes of operation, each corresponding to a system of LDIs. Let us take the example of the flow shop again. By assigning a mode of operation to each job type, we can model the manufacturing system by SLDIs in which different job entrance orders simply correspond to different schedules, i.e., sequences of modes. A first advantage of SLDIs is thus that, with a single dynamical system, one can represent all possible entrance orders of jobs in the flow shop, including non-periodic ones.
It is worth noting that this property is not exclusive to SLDIs, as also Ptime Petri nets can be utilized to model flow shops under any job order [19,20]. In fact, there exists a strong relation between the two models (which is briefly discussed in Section 6). However, in contrast to P-time Petri nets, the dynamics of SLDIs possesses another appealing feature: a switched-linear formulation in the max-plus algebra. By exploiting the abundance of existing results in this algebraic framework (see, e.g., [21][22][23][24]), in this paper we will derive lowcomplexity algorithms for the cycle time computation, considering two types of schedules: periodic and intermittently periodic schedules. In the example of the flow shop, the first type of schedules corresponds to periodic arrivals of jobs of different types in the manufacturing system (as in schedule ababab . . .), whereas the second one is a generalization in which the entrance order of jobs alternates among periodic and non-periodic regimes. An example of a schedule of the second type for the flow shop of Figure 1 is aababab . . . abbbababa . . ., which consists of an initial transient a followed by periodic regime abab . . ., an intermediate transient b, and the final periodic subschedule baba . . . 3 In the cycle time analysis for this kind of schedules, one is interested in finding the cycle times that can be achieved in each periodic regime. The motivation for studying intermittently periodic schedules is two-fold: 1. they are ubiquitous in applications (as discussed in Section 4.4), and 2. differently from systems described by only lower time-bounds, the cycle time analysis in this class of schedules can produce non-trivial results in the presence of time-window constraints. Here, by "non-trivial" we mean that the cycle times of the periodic subschedules in an intermittently periodic schedule may be different from those obtained by studying the periodic subschedules independently.
The present article enhances and extends the recent conference paper [25] in several ways. Besides simplifying the preliminaries in Section 2 and adding a number of simple and illustrative examples, an important contribution of this extended paper is to systematically investigate the initial conditions of P-TEGs (in Section 3) and their relation to SLDIs. Two types of initial conditions are presented, loose and strict, and it is proven that P-TEGs with strict initial conditions can be represented by SLDIs but not by pure LDIs (Section 4.2). In Section 4, after formally presenting SLDIs, the cycle time analysis in periodic and intermittently periodic trajectories is discussed; intermittently periodic trajectories are introduced for the first time in this extended version. In [25], the correctness of the low-complexity method for computing the cycle time in SLDIs under periodic schedules was proven using tools from automata and regular languages theory; in this paper, the proof is entirely based on algebraic arguments. Although the initial inspiration for the algorithm was drawn from analogies between multi-precedence graphs and automata theory, the new proof relies on more established results, which hopefully makes it less arduous to read.
In Section 5, the analysis of the case study considered in [25] has been deepened further. The case study consists of a robotic job shop example derived from [26], where parts of different type are required to visit a sequence of processing stations in different order, and are transported by a single robot. The authors of [26] proved that the cycle time analysis in this class of systems can be performed in strongly polynomial time complexity O(V 4 n 4 ), where V is the period of the entrance order of different types of parts in the system, and n is the number of processing stations. In this paper, we show that the complexity can be reduced to O(V n 3 + n 4 ) using SLDIs. Computational tests show that the advantage is not only theoretical, but translates into a tangibly faster cycle time analysis. This makes our algorithm an appealing optimization subroutine for the solution of (NP-hard) cyclic scheduling problems in which the goal is to find the optimal path of the robot. Additionally, in this paper we show how to derive a complete (intermittently periodic) trajectory of the system, consisting of a start-up transient (where parts are initially introduced into the system), a periodic regime, and a shut-down transient (in which all parts are removed from the stations).
Finally, Section 6 provides concluding remarks, comparisons with related classes of dynamical systems, and suggestions for future work.

Notation
The set of positive, respectively non-negative, integers is denoted by N, respectively N 0 . The set of non-negative real numbers is denoted by R ≥0 . Moreover, R max := R ∪ {−∞}, R min := R ∪ {∞}, and R : n×n , we will use notation A ♯ to indicate −A ⊤ . Given a, b ∈ Z with b ≥ a, a, b denotes the discrete interval {a, a + 1, a + 2, . . . , b}.

Preliminaries
In this section, some basic concepts of max-plus algebra (Section 2.1) and precedence graphs (Section 2.2) are recalled; for a more detailed discussion of those topics, we refer to [22][23][24]. Sections 2.3 and 2.4 present the nonpositive circuit weight problem (introduced in [27]) and max-plus linear-dual inequalities.

Max-plus algebra
The max-plus algebra operates on the set of real numbers extended with −∞ and +∞, and is endowed with operations ⊕ (addition), ⊗ (multiplication), ⊞ (dual addition), and ⊠ (dual multiplication), defined as follows: for all a, b ∈ R, These operations can be extended to matrices; given A, B ∈ R m×n , C ∈ R n×p , When the meaning is clear from the context, we will denote the max-plus multiplication between matrices A and C, A ⊗ C, simply by AC. The symbols E , T , and E ⊗ denote, respectively, the neutral element for ⊕, ⊞, and ⊗, i.e., E ij = −∞ and T ij = +∞ for all i, j, and E ⊗ is a square matrix with ( Given a square matrix A, its r th power is defined recursively by A 0 = E ⊗ and, for all r ≥ 1, A r = A r−1 ⊗ A. Moreover, the Kleene star of a matrix A ∈ R n×n is The partial order relation ≤ between two matrices A and B of the same size is defined elementwise: A ≤ B if and only if A ij ≤ B ij for all i, j. Analogously to the standard algebra, we define the max-plus trace of matrix A ∈ R n×n by The product and dual product between scalar λ ∈ R and matrix A ∈ R m×n are given by If λ / ∈ {−∞, +∞}, the two expressions coincide. We will therefore simply write λA in place of λ⊗A or λ ⊠ A when λ ∈ R. When λ ∈ R, we indicate by λ −1 the element such that λ −1 ⊗ λ = λ ⊗ λ −1 = 0; thus, in the standard algebra, λ −1 coincides with −λ. When λ ∈ {−∞, +∞}, λ does not have a multiplicative inverse; nevertheless, we will use symbol λ −1 to denote −λ.

Precedence graphs
A directed graph is a pair (N, E) where N is a finite set of nodes and E ⊆ N ×N is the set of arcs. A weighted directed graph is a triplet (N, E, w), where (N, E) is a directed graph, and w : E → R is a function that associates a weight w((i, j)) to each arc (i, j) ∈ E of graph (N, E). A sequence of r + 1 nodes ρ = (i 1 , i 2 , . . . , i r+1 ), r ≥ 0, such that (i j , i j+1 ) ∈ E for all j ∈ 1, r is a path of length r; a path ρ such that i 1 = i r+1 is called a circuit. The weight of a path is the sum (in conventional algebra) of the weights of the arcs composing it; conventionally, the weight of a path of length r = 0, i.e., ρ = (i 1 ), is equal to 0.
The precedence graph associated with matrix A ∈ R n×n max is the weighted directed graph G(A) = (N, E, w), where N = 1, n , and there is an arc (j, i) ∈ E of weight w((j, i)) = A ij if and only if A ij = −∞. We say that G(A) is a parametric precedence graph when elements of A are functions of some real parameters λ 1 , . . . , λ p , i.e., A = A(λ 1 , . . . , λ p ).
There are important connections between the max-plus algebra and precedence graphs. For instance, element (i, j) of the r th max-plus power of a matrix A ∈ R n×n max , (A r ) ij , corresponds to the maximum weight of all paths in G(A) of length r from node j to node i. A direct consequence is that (A * ) ij is equal to the largest weight of all paths (of any length) from node j to i. Observe that a precedence graph G(A) does not contain circuits with positive weight if and only if tr(A * ) = 0; in presence of at least one circuit with positive weight, tr(A * ) = ∞. In the following, we indicate by Γ the set of all precedence graphs that do not contain circuits with positive weight, i.e., The following proposition will be used later to verify the existence of, and compute, trajectories satisfying time-window constraints in (switched) maxplus linear-dual inequalities.
Proposition 1 [22,23] Let A be an n × n matrix with elements in Rmax. Inequality A ⊗ x ≤ x admits a solution x ∈ R n if and only if G(A) ∈ Γ. In this case, any column of A * solves the inequality, i.e., A ⊗ (A * ) ·,i ≤ (A * ) ·,i for all i ∈ 1, n .
The maximum circuit mean of precedence graph G(A), denoted by mcm(A), is the maximum weight-over-length ratio of all circuits of positive length in the graph; this value coincides with the largest max-plus eigenvalue (the max-plus spectral radius) of A, i.e., the largest λ ∈ R max such that A ⊗x = λx for some vector x with elements from R max . For a matrix of dimension n × n, the maximum circuit mean can be computed through the following formula ( [22]): where a 1 k (corresponding to a k in the standard algebra) is the k th max-plus root of a ∈ R max ; a more efficient algorithm that returns the same value in time complexity O(n × m) in the worst case, where m is the number of edges in G(A), is due to Karp [28].
We recall that it is possible to check whether G(A) ∈ Γ and, when G(A) ∈ Γ, to compute A * in time O(n 3 ) in the worst case, using the Floyd-Warshall algorithm [29]. In the case G(A) / ∈ Γ, computing A * is an NP-hard problem; fortunately, the algorithms presented in the next sections will never face this issue in practice.

The non-positive circuit weight problem
Given a parametric precedence graph G(A), where A = A(λ 1 , . . . , λ p ), the non-positive circuit weight problem (NCP) consists in characterizing the set Λ NCP (A) = {(λ 1 , . . . , λ p ) ∈ R p | G(A) ∈ Γ} of all values of parameters λ 1 , . . . , λ p for which G(A) does not contain circuits with positive weight. Specific classes of the NCP find applications in the analysis of periodic trajectories in max-plus dynamical systems. An application example is shown in Section 2.4.

The PIC-NCP
When matrix A has the form A(λ) = λP ⊕ λ −1 I ⊕ C for arbitrary matrices P, I, C ∈ R n×n max (called proportional, inverse 4 , and constant matrix, respectively), then the problem is referred to as the proportionalinverse-constant-NCP. In this case, Λ NCP (λP ⊕ λ −1 I ⊕ C) = [λ min , λ max ] ∩ R is an interval, and its extreme points can be found either in weakly polynomial time using linear programming solvers such as the interior-point method, or in strongly polynomial time O(n 4 ) using Algorithm 1 ( [27]). The functioning of the latter algorithm is briefly described in the following.
We remark that the PIC-NCP represents a generalization of the max-plus subeigenproblem (see [30]), i.e., the problem of finding a real λ such that the inequality I ⊗ x ≤ λx admits a solution x ∈ R n for a given matrix I ∈ R n×n max . Indeed, the above inequality can be rewritten (by multiplying both sides by λ −1 ) as and, from Proposition 1, admits a solution x ∈ R n if and only if the precedence graph G(λ −1 I) does not have circuits with positive weight; the PIC-NCP thus simplifies into the max-plus subeigenproblem when matrices P and C are E. We recall from [30, Lemma 1] that the least solution of the subeigenproblem is the max-plus spectral radius of matrix I, i.e., Λ NCP (λ −1 I) = [mcm(I), +∞) ∩ R.
When matrices P and C are not E, a more sophisticated approach is necessary to solve the PIC-NCP. In particular, in Algorithm 1 some pre-computations are first performed (lines 3-4) to simplify the problem into an equivalent PI-NCP (a PIC-NCP where matrix C is E) by redefining P as C * ⊗ P ⊗ C * and I as C * ⊗ I ⊗ C * . Then, through the for-loop of lines 5-7, a matrix S is constructed such that, in every new iteration of the loop, the spectral radius of matrix I ⊗ S * and the inverse of the spectral radius of matrix P ⊗ S * approximate better and better λ min and λ max , respectively. It can be shown (see [27]) that after at most n 2 iterations, the two quantities converge to the desired values (line 10), i.e., Algorithm 1: Solve NCP(P, I, C) (from [27]) If some conditions on matrices C and S do not hold, the algorithm can be terminated prematurely as no λ solving the PIC-NCP exists (lines 1-2 and 8-9).

The MPIC-NCP
A natural generalization of the PIC-NCP is the multivariable PIC-NCP (MPIC-NCP), in which matrix A takes the form for some matrices P i , I i , C ∈ R n×n max and parameters λ i ∈ R, for all i ∈ 1, q . In this case, it can be shown that the problem becomes "trivially intractable" in q, 5 in the sense that the set Λ NCP (A) corresponds, in the worst case, to the solution of a system of 3 q − 1 (i.e., exponentially many) non-redundant linear inequalities (in the conventional sense) in variables λ 1 , . . . , λ q ; in other words, Λ NCP (A) is a polytope with at most 3 q − 1 facets living in a q-dimensional space. Consequently, it is unrealistic to expect to efficiently solve the MPIC-NCP, since the solution set Λ NCP (A) cannot be described in polynomial space. Nevertheless, it is possible to verify the non-emptiness of Λ NCP (A) in (weakly) polynomial time from the following observation. Proposition 1 suggests that the parameters λ i such that G(A) ∈ Γ are those for which the inequality admits a solution x ∈ R n . We can rewrite the inequality above in the standard algebra as the system which is equivalent to (1) The system above consists of at most 6 (2q+1)n 2 linear inequalities in n+q real unknowns x 1 , . . . , x n , λ 1 , . . . , λ q . Therefore, the non-emptiness of its solution set can be checked in polynomial time using linear programming techniques.

Max-plus linear-dual inequalities
In the following, we define max-plus linear-dual inequalities (LDIs) from a purely formal point of view. Their application to describe the dynamics of P-time event graphs is discussed in the next section. Let A 0 , A 1 ∈ R n×n max , B 0 , B 1 ∈ R n×n min , and K ∈ N ∪ {+∞}. LDIs are systems of (⊕, ⊗)and (⊞, ⊠)linear dynamical inequalities in the dater function x : 1, K → R n of the form ∀k ∈ 1, K , A finite (when K ∈ N) or infinite (when K = +∞) trajectory {x(k)} k∈ 1,K of length K is consistent if it satisfies (2) for all k. It is often useful in practice to restrict the focus on the simple class of 1-periodic trajectories, which are those of the form {λ k−1 x(1)} k∈ 1,K ; in the standard algebra, 1-periodic trajectories are those that satisfy: ∀k ∈ 1, K , i ∈ 1, n , x i (k) = (k − 1) × λ + x i (1). The number λ is called period or cycle time of the 1-periodic trajectory.
In the following, we recall how to verify the existence of 1-periodic trajectories for given LDIs. Substituting in (2) the formula x(k) = λ k−1 x(1), we obtain ∀k ∈ 1, K , which, after multiplying left and right hand sides of the above inequalities by (λ k−1 ) −1 , 7 simplifies to We recall the following result.
Because of Proposition 2 we can rewrite (3) as Note that we obtained a PIC-NCP with matrices P = B 1♯ , I = A 1 , and C = A 0 ⊕ B 0♯ . Therefore, a consistent 1-periodic trajectory exists if and only if Λ where λ min and λ max can be found in time O(n 4 ) using Algorithm 1.

P-time event graphs
Definition 1 (From [33]) An ordinary (or unweighted) P-time Petri net (P-TPN) is a 5-tuple (P, T , E, m, ι), where (P ∪ T , E) is a directed graph in which the set of nodes is partitioned into the set of places, P, and the set of transitions, T , the set of arcs, E, is such that E ⊆ (P × T ) ∪ (T × P), m : P → N 0 is a map such that m(p) represents the number of tokens initially residing in place p ∈ P (also called initial marking of p), and ι : In the following, we briefly describe the dynamics of an ordinary P-TPN 9 . A transition t is enabled when either it has no upstream places or each upstream 7 Recall that, in the standard algebra, this is equivalent to subtracting everywhere by (k−1)×λ. 8 The reader familiar with residuation theory [31] may note that function f ♯ A (y) = A ♯ ⊠ y coincides with the residual of the function "left multiplication by A" fA(x) = A ⊗ x (i.e., f ♯ A (y) is the greatest x that satisfies A ⊗ x ≤ y). In the max-plus context, function f ♯ A is often called "left division by A" and denoted by f ♯ A (y) = A • \y [22]. Dually, function g ♭ B (x) = B ♯ ⊗ x coincides with the dual residual of the function "left dual multiplication by B" gB(y) = B ⊠ y (i.e., g ♭ B (x) is the least y that satisfies x ≤ B ⊠ y), and is sometimes denoted by g ♭ B (x) = B • \x [32]. 9 Since the functioning of P-TPNs is here reported only to present P-time event graphs, we will not delve into the nuances of weak and strong semantics for P-TPNs. For a thorough description of different types of semantics, we invite the reader to consult [34]. place p of t contains at least one token which has resided in p for a time between τ − p and τ + p (extremes included). When transition t is enabled, it may fire; its firing causes one token to be removed instantaneously from each of the upstream places of t, and one token to be added, again instantaneously, to each of the downstream places of t. If a token sojourns more than τ + p time instants in a place p, then the token becomes dead.
A P-time event graph (P-TEG) is an ordinary P-TPN in which every place has exactly one upstream and one downstream transition. Let |T | = n be the number of transitions in a P-TEG and let x : 1, K → R n be a dater function of length K ∈ N ∪ {+∞}, i.e., a function such that x i (k) represents the time at which transition t i fires for the k th time. Since the (k + 1) st firing of any transition cannot occur before the k th , we require the dater to be a non-decreasing function, i.e., ∀i ∈ 1, n , x i (k + 1) ≥ x i (k). The evolution of the marking in a P-TEG is entirely described by its corresponding dater trajectory {x(k)} k∈ 1,K , and if a non-decreasing dater trajectory exists for which no token death occurs, then it is said to be consistent for the P-TEG.
It is always possible to transform a P-TEG into an equivalent one whose places have at most 1 initial token each [17]. Therefore, in the following we will only focus on P-TEGs in which the initial marking m(p) is either 0 or 1 for each place p ∈ P. Under this assumption, a consistent trajectory for a given P-TEG must satisfy LDIs as in (2), where matrices A 0 , A 1 ∈ R n×n max , B 0 , B 1 ∈ R n×n min are called characteristic matrices of the P-TEG, and are defined as follows. If there exists a place p ij with initial marking µ ∈ {0, 1}, upstream transition t j and downstream transition t i , then A µ ij = τ − pij and B µ ij = τ + pij ; otherwise, A µ ij = −∞ and B µ ij = ∞. Before introducing some structural properties of P-TEGs, it is useful to clarify the role of initial conditions for these dynamical systems.

Initial conditions
Depending on the P-TEG's intended application, the restrictiveness of initial conditions may vary; in the following, we introduce two alternatives.

Loose initial conditions
Inequalities (2) suggest that x(1) can assume any value in R n , as long as it satisfies Note that (4) does not restrict the first firing time based on the arrival times of the initial tokens in the Petri net; the arrival times of initial tokens are, in fact, not even defined. For this reason, we say that the initial conditions of a P-TEG are loose if no other restriction other than (4) applies to x(1). P-TEGs with loose initial conditions evolve entirely according to LDIs, and are suitable to model manufacturing systems operating in periodic regime, after a transient period has passed (an example is given in Section 5.1). Another   scenario where loose initial conditions may be convenient is when the timewindow constraints need to be fulfilled only after the occurrence of the first events, as shown in the following example.
Example 1 (Heat treatment line) This example is adapted from [35]. Consider a simple heat treatment line, schematically represented in Figure 2, consisting of a furnace, which performs a heat treatment, and an autonomous guided vehicle that receives processed pieces and transports them to the next stage. Both the furnace and the vehicle have unitary capacity, i.e., they can process and transport one piece at a time, respectively. The heat treatment must last between 2 and 3 time units, and the autonomous vehicles takes (at least) 0.5 time units both to transport a processed piece to the next stage and to travel back to the furnace. Customers' demand imposes that the time difference between subsequent unloadings of processed pieces from the autonomous guided vehicle must not exceed 4 time units; the specification needs to be met for all pieces after the first one. Moreover, each piece must spend at least 6 time units in the processing line, from the moment it enters the furnace to the one it is removed from the vehicle, in order to synchronize with other processing stages. Initially, the furnace is empty and the vehicle is waiting for a piece at the furnace. The P-TEG in Figure 3 models the described plant; a firing of transitions t 1 , t 2 , and t 3 represents, respectively, the arrival of an unprocessed piece in the furnace, the loading of a processed piece onto the autonomous guided vehicle, and the unloading of a piece from the vehicle.
The characteristic matrices of the P-TEG are It is possible to verify that is a consistent trajectory for the P-TEG under loose initial conditions. Observe that the first firing time of transition t 3 does not violate the upper bound associated to place p 33 (i.e., the place that is upstream and downstream of transition t 3 ), even though x 3 (1) = 6 > 4 = B 1 33 . Indeed, the sojourn time of the initial token in place p 33 does not restrict the dynamics of the P-TEG. This is convenient from a practical point of view, as the constraint on the processing rate must be enforced only after the first piece leaves the plant.

Strict initial conditions
For the considered application, it may be necessary to impose further restrictions on the initial conditions. Here we take into account the amount of time that initial tokens have resided in places prior to the initial time t 0 ∈ R; we call this value the time tag of the token 10 . Time tags can be useful, for instance, in manufacturing, to specify that some machines have already been processing a part since time t 0 − τ , or in transportation, to indicate that a vehicle has left a station at time t 0 − τ , for τ ≥ 0.
Let ρ be a function that associates a time tag to every place with an initial token in a P-TEG. Formally, if there is a place p ij with marking m(p ij ) = 1, upstream transition t j , and downstream transition t i , then we denote ρ(p ij ) = ρ ij ∈ R ≥0 , otherwise ρ(p ij ) is not defined. Then, in addition to (2), the first firing time of the transitions of the P-TEG must satisfy, for all i, j ∈ 1, n , for each i, j, the inequality specifies that the first firing time of transition t i must not violate the time-window constraint [A 1 ij , B 1 ij ] associated with place p ij , considering that the initial token of this place arrived at time t 0 − ρ ij . In the max-plus algebra, the latter inequalities can be expressed as  where Note that other possible definitions for ∆ and ∆, corresponding to different requirements for the first firings of transitions, may be considered. In general, given any ∆ ∈ R n×n max , ∆ ∈ R n×n min such that (∆, ∆) = (E , T ), inequality (5) restricts the set of consistent trajectories for a P-TEG; hence, we say that the initial conditions of a P-TEG are strict if x(1) is required to satisfy them for some (∆, ∆) = (E , T ). We will refer to consistent trajectories with either loose or strict initial conditions depending if they satisfy only (2) or also (5). Note that, without loss of generality, we can assume that t 0 = 0, as P-TEGs (with either loose or strict initial conditions) are time-invariant systems, i.e., if {x(k)} k∈ 1,K is a consistent trajectory, then {t 0 ⊗ x(k)} k∈ 1,K is consistent as well for any t 0 ∈ R. In other words, the choice of the initial time t 0 does not affect the dynamics of P-TEGs.
Example 2 (Heat treatment line, cont.) Consider again the P-TEG of Figure 3. It is not difficult to see that, if we assign a time tag to each initial token of the P-TEG, no consistent trajectory that satisfies strict initial conditions can be found; indeed, it is not possible to fire transition t 3 before time 4 − ρ 33 , for any time tag ρ 33 ∈ R ≥0 . So, let us modify the configuration of the initial tokens as in Figure 4; time tags are indicated in the figure.
The interpretation is that: • a piece is inside the heat treatment line since time t 0 − 3, as ρ 31 = 3, • an autonomous guided vehicle is at the unloading location with a processed piece from time t 0 , as ρ 32 = 0.5 = A 1 32 , • the furnace has completed the last heat treatment at time t 0 −0.5, as ρ 12 = 0.5, and • the first processed piece is required to leave the heat treatment plant before time t 0 + 3, as ρ 33 = 1.
The characteristic matrices for this example are and the matrices ∆ and ∆ are Assuming that t 0 = 0, the following is a consistent trajectory for the P-TEG under strict initial conditions: Despite their usefulness in applications, strict initial conditions present an additional complexity: the inequalities describing the dynamics of of P-TEGs with strict initial conditions are not (pure) LDIs. This means that mathematical results in LDIs can be directly applied to P-TEGs with loose but not with strict initial conditions; this will be made evident in the following section. In Section 4.2, we will see that the dynamics of P-TEGs with strict initial conditions falls in the category of switched LDIs.

Structural properties
In this section we recall the definition of some structural properties of P-TEGs. These properties can be equivalently stated for P-TEGs under loose or strict initial conditions.
A P-TEG is said to be consistent if it admits a consistent, non-decreasing trajectory {x(k)} k∈N of infinite length. The non-decreasingness of the dater trajectory, equivalent to having x(k + 1) ≥ x(k) for all k ∈ 1, K − 1 , is a natural requirement, as the (k + 1) st firing of a transition cannot occur before the k th one; because of Proposition 2, to restrict the evolution of the dater trajectory such that this restriction is always satisfied, it is sufficient to modify the definition of matrix  Figure 5: Example of boundedly consistent P-TEG with strict initial conditions that admits no 1-periodic trajectory.
We say that a trajectory {x(k)} k∈N is delay-bounded if there exists a positive real number M such that, for all i, j ∈ 1, n and for all k ∈ N, x i (k) − x j (k) < M ; a P-TEG admitting a consistent delay-bounded trajectory of the dater function is said to be boundedly consistent. Although in consistent P-TEGs it is possible to find a marking evolution such that no time-window constraint is violated, if the stronger property of bounded consistency does not hold, any consistent, infinite trajectory will accumulate unbounded delay between the firing times of two distinct transitions. This phenomenon is usually not desirable in manufacturing systems represented by P-TEG, where the firings of transitions represent the start or end of processes, and the k th product entering the system is finished when all transitions fire for the k th time. Indeed, in this context it implies that the total time the k th product spends in the manufacturing system increases without bounds with k.
Analogously to LDIs, in P-TEGs we say that dater trajectories of the form {λ k−1 x(1)} k∈ 1,K are 1-periodic with period λ ∈ R ≥0 . Clearly, in P-TEGs with loose initial conditions, 1-periodic trajectories can be found in time complexity O(n 4 ) using Algorithm 1, as their evolution satisfies LDIs. To our knowledge, no algorithm that checks whether a P-TEG is consistent has been found until now; on the other hand, bounded consistency of P-TEGs with loose initial conditions can be verified in time O(n 4 ). This fact comes from the following result.

Theorem 3 ([35]) A P-TEG with loose initial conditions is boundedly consistent if and only if it admits a consistent 1-periodic trajectory.
On the other hand, an analogous result for the case with strict initial conditions does not hold: boundedly consistent P-TEGs with strict initial conditions may admit no 1-periodic trajectory, as shown in the following example.
Example 3 Using an algorithm that will be presented in Section 4.2, it can be shown that the P-TEG with strict initial conditions in Figure 5 admits no 1-periodic trajectory. However, assuming that t 0 = 0, it admits the following delay-bounded (2-periodic) trajectory: Therefore, it is boundedly consistent.
The following example illustrates the discussed properties in the case of P-TEGs with loose initial conditions. Example 4 Consider the P-TEG represented in Figure 6, in which time windows are parametrized with respect to label z; in Table 1, values of time windows are given for z ∈ {a, b, c}. The matrices characterizing the P-TEG labeled z are: We analyze structural properties of the P-TEGs under loose initial conditions.  Table 1: Parameters for the P-TEG of Figure 6.
Since lower and upper bounds for the sojourn times of the two places with an initial token coincide, once the vector of first firing times xz(1) is chosen (such that the first inequality in (2) is satisfied for k = 1, i.e., x z,2 (1) ≥ x z,1 (1)), the only infinite trajectory {xz(k)} k∈N that is a candidate to be consistent for the P-TEG labeled z is deterministically given by However, for the case z = a it is easy to see that, for any valid choice of the vector of first firing times, the candidate trajectory {xa(k)} k∈N is not consistent (as for a sufficiently large k, x a,2 (k) < x a,1 (k), i.e., the first inequality of (2) is violated). For z = b, candidate trajectories {x b (k)} k∈N , despite being consistent, are not delaybounded and result in the infinite accumulation of tokens in the place between t 1 and t 2 for k → ∞. On the other hand, {xc(k)} k∈N is consistent and delay-bounded (in fact, it is 1-periodic with period 1). Thus we can conclude that the P-TEG labeled a is not consistent, the one labeled b is consistent but not boundedly consistent, and the one labeled c is boundedly consistent. Of course, we would have reached the same conclusion regarding bounded consistency by using Theorem 3.

Mathematical description
We start by defining switched LDIs (SLDIs) as the natural extension of LDIs in which matrices A 0 , A 1 , B 0 , B 1 may be different for all k. Formally, SLDIs are a 5-tuple S = (Σ, A 0 , A 1 , B 0 , B 1 ), where Σ = {a 1 , . . . , a m } is a finite alphabet whose symbols are called modes, and A 0 , A 1 : Σ → R n×n max , B 0 , B 1 : Σ → R n×n min are functions that associate a matrix to each mode of Σ; for the sake of simplicity, given a mode z ∈ Σ, we will write A 0 , respectively. We denote by Σ * and Σ ω the sets of finite and infinite concatenations of modes from Σ, respectively. A schedule w is an element of Σ * ∪ Σ ω , i.e., it is either a finite or an infinite sequence of The dynamics of SLDIs S under schedule w ∈ Σ * ∪ Σ ω is expressed by the following system of inequalities: where function x : 1, K → R n is called dater of S associated with schedule w. Term x i (k) represents the occurrence time of event i associated with mode 11 w k . Similar to P-TEGs, it is natural to assume the following nondecreasingness condition for the dater of SLDIs: for all k, h ∈ 1, K , h > k, such that w k = w h , x(h) ≥ x(k). The implication is that events occurring during a later execution of a mode cannot occur before events that took place during an earlier execution of that mode. Note that this does not imply that x(k) is non-decreasing over k; this stronger condition would indeed unnecessarily limit the modeling expressiveness of SLDIs, as illustrated by the dater trajectory in Example 6 at page 21.
For convenience, given a finite sequence of modes v = v 1 v 2 . . . v V ∈ Σ * of length V ∈ N and a number K ∈ N, in the remainder of the paper we will denote by v K ∈ Σ * the sequence of length V · K formed by concatenating We now show possible applications of SLDIs with two simple examples. were to be processed. The case where only parts of type a are to be processed is shown in Figure 3.
Example 5 (Heat treatment line, cont.) Consider again the heat treatment line of Example 1. Now, suppose that two types of parts can be processed in the system: part a and part b; in this example, a schedule w ∈ Σ * ∪ Σ ω represents the entrance order of parts in the heat treatment line. As illustrated in Figure 7, the two parts require different heating times; pieces of type a must be heated in the furnace for a time between 2 and 3 time units (as in Example 1), whereas pieces of type b between 3 and 4 time units. Moreover, the processing rate requirement changes depending on the type w k ∈ {a, b} of the k th part entering the plant: the (k + 1) st part must be unloaded from the autonomous guided vehicle at most 4 time units after the k th one if w k = a, and at most 5 time units later if w k = b. As in Example 1, we consider loose initial conditions, i.e., we suppose that the processing rate requirement needs to hold for all pieces after the first one. This system can be modeled by SLDIs Clearly, if w = aaa . . . = a K is a concatenation of only mode a, the dynamics of the system can be described by the P-TEG of Figure 3; similarly, if w = b K , then the SLDIs are equivalent to the dynamics of the P-TEG of Figure 8.
In this simple example, the following trajectory satisfies all the time-window constraints, for any schedule w ∈ {a, b} * ∪ {a, b} ω : The example shows that SLDIs are capable of modeling flow shops with time-window constraints, i.e., manufacturing systems where different jobs (in this case, parts) are processed in each machine (the furnace and the autonomous guided vehicle) of the system in the same order.
Example 6 (Starving philosophers problem) We present a variant of the famous dining philosophers problem, which we call the starving philosophers problem. There are p ∈ N philosophers sitting at a table eating spaghetti; on the table, there are p chopsticks, and each philosopher i ∈ 1, p − 1 needs both the i th and the (i + 1) st chopstick to start eating, whereas the p th philosopher needs the p th and the 1 st chopstick. The i th philosopher takes c ij time units to grab the j th chopstick and e i time units to eat; philosophers are allowed to grab two chopsticks at the same time and are not forced to stop eating after e i time units. After eating, a philosopher instantaneously puts the chopsticks on the table, so that they can be used by other philosophers. In our version of the problem, dining can take a "macabre" turn: if the i th philosopher does not eat for more than s i time units, s/he will starve. The objective of the problem is to find a dining order such that no philosopher starves. The problem is a metaphor for an issue encountered in concurrent programming, namely, resource starvation. Consider p critical processes (the philosophers) that need to access some shared resources (the chopsticks); for safety reasons, it might be desirable to prevent some processes from not receiving the requested resources for too long. Thus, a scheduling plan should guarantee safe operation by granting each process the permission to access the resources at the right time.
Here we suppose that there are p = 4 philosophers at the table, and that the following periodic dining order is imposed: initially, the second and the fourth philosophers eat (they can do so concurrently, as they do not need to share chopsticks); after that, the first philosopher eats once while, in the meantime, the third philosopher eats twice in a row; finally, the eating order repeats from the beginning. The order in which philosophers eat before repeating the periodic sequence is referred to as the dining cycle.
Since the considered dining order is periodic, it is possible to describe all trajectories that are valid for the system by means of a P-TEG (with time tags, if we suppose that the i th philosopher should start eating for the first time before time t 0 + s i for all i ∈ 1, p ). The P-TEG for this example is shown in Figure 9, where a firing of transitions t s,i and t f,i represents that the i th philosopher has started and finished eating for the first time in a dining cycle, and a firing of transitions t ′ s, 3 and t ′ f,3 indicate that the third philosopher has started and finished eating for the second time in a dining cycle, respectively. Note that the dimension of the dater function for this problem increases not only with the number of philosophers, but also with the amount of times a philosopher eats in a dining cycle; furthermore, observe that P-TEGs can represent infinite eating orders only if they are periodic.
The same system can be represented more compactly by SLDIs. Let Σ = {i, p 1 , p 2 , p 3 , p 4 } be the alphabet associated with the SLDIs, where i is the auxiliary initial mode, which will be used to impose strict initial conditions on the system (strict initial conditions are analyzed in more depth in Section 4.2), and p i corresponds to the i th philosopher. A meaningful schedule for this system is any sequence w ∈ Σ K such that w 1 = i and, for all k ∈ 2, K , w k = p i k for some i k ∈ 1, p . Whereas the first mode w 1 = i does not have physical interpretation besides mathematically characterizing the initial conditions for the system, for all k ∈ 2, K , w k describes the eating order of philosophers. The interpretation of schedule w is as follows; consider k ∈ 2, K , w k = p i and w k+1 = p j : • if the i th and j th philosophers require access to the same chopstick, then philosopher i will eat once before philosopher j; • otherwise, the i th and j th philosophers will eat independently of each other, i.e., philosopher i will start (and finish) eating either before or after or at the same time as philosopher j. In this case schedules uw k w k+1 v and uw k+1 w k v are representative of the same behavior of the system 12 , for u ∈ Σ * and v ∈ Σ * ∪Σ ω .
A possible schedule corresponding to the chosen dining order is then given by We will design A 0 , A 1 , B 0 , B 1 such that any dater function x(k) ∈ R p+1 = R 5 satisfying (6) assumes the following interpretation, from which the evolution of the system can be obtained: for all k ∈ 2, K , if w k = p i , then x i (k) and x 5 (k) represent, respectively, the time at which the i th philosopher starts and finishes eating; therefore, assuming w k+1 = p j , if both the i th and the j th philosophers require access to the h th chopstick, then x i (k) + e i + c jh ≤ x j (k + 1) (sequential behavior), otherwise, if they do not need to share chopsticks, x i (k) could also be greater than x j (k + 1) (concurrent behavior). For all philosophers i ∈ 1, 4 such that w k = p i , x i (k) is an auxiliary variable that stores the time at which the i th philosopher will eat next. For all i ∈ 1, 5 , element x i (1) will be assigned to the initial time t 0 , in order to manage the initial conditions (for more details, see Section 4.2). For example, consider a value of k ∈ 2, K − 1 such that w k = p 1 ; with the above interpretation, in order to represent the dynamics of the system, the dater function must satisfy where (7a) imposes the time between starting and finishing eating for the 1 st philosopher, (7b) is used to force the 1 st philosopher to start eating again only after s/he has grabbed once more the 1 st and 2 nd chopsticks but before starving, (7c) and (7d) impose the 2 nd and 4 th philosophers to start eating only after grabbing the chopsticks left by the 1 st philosopher, and (7e-7g) are auxiliary constraints to impose that x i (k + 1) = x i (k) for all philosophers i ∈ 2, 4 . Finally, for the initial condition, we want to impose that max(c 11 , c 12 ) (2) ≤ s 4 + t 0 to make sure that philosophers start eating for the first time after grabbing the chopsticks and before starving.
The first 5 elements of the dater trajectory for the SLDIs are: It is worth noting that, differently from P-TEGs, the dater function of SLDIs does not need to be non-decreasing: for instance, in this example we have x 5 (4) = 7.5 ≥ 6 = x 5 (5), as the time in which the first philosopher (w 4 = p 1 ) stops eating for the first time is after the time in which the third philosopher (w 5 = p 3 ) stops eating for the first time.
The three main advantages of using SLDIs instead of P-TEGs for this problem are the following: [c23, ∞] Figure 9: P-TEG for the starving philosophers problem. Tokens inside a place colored black, blue, and red represent, respectively, a philosopher eating, a philosopher waiting to eat, and a chopstick being grabbed. The time tag associated to each place with initial token is 0. Note that a token in the place with upstream transition t f,3 and downstream transition t ′ s,3 indicates that the third philosopher is both grabbing the third and the fourth chopstick and waiting to eat, hence the double color of this place.
1. higher computational efficiency: the dater function for the SLDIs has smaller dimension compared to the P-TEG. As we shall see, this corresponds to lower computational complexity for analyzing trajectories of the system; 2. lower modeling effort: the P-TEG in Figure 9 can only represent the dining order specified above; to analyze a different dining order, a new P-TEG needs to be provided. On the other hand, for the SLDIs different dining orders simply correspond to different schedules w; 3. larger modeling expressiveness: only SLDIs are able to represent dining orders that are not periodic (with a dater function of finite dimension).
When schedule w is fixed, we can extend the definition of some properties of P-TEGs to SLDIs in a natural way. For instance, if there exists a trajectory of the dater {x(k)} k∈ 1,K that satisfies (6), then the trajectory is consistent for the SLDIs under schedule w, and we say that w is a consistent schedule for the SLDIs (or that the SLDIs are consistent under schedule w). Note that, different from the simple case of Example 5, there are SLDIs for which not all schedules admit consistent trajectories.
The definitions of delay-bounded trajectory and bounded consistency are generalized to SLDIs in a similar fashion. The interpretation of bounded consistency of a schedule w is analogous to the one of P-TEGs; when a process consisting of several tasks (the start and end of which are represented by events) is modeled by SLDIs under a schedule w that is not boundedly consistent, then either the execution of every possible sequence of tasks following w will lead to the violation of some time window constraints (if w is not even consistent), or we will certainly observe an infinite accumulation of delay between the start or end of some tasks (if the only consistent trajectories are not delay-bounded).

SLDIs and P-TEGs with strict initial conditions
As discussed in Subsection 3.1.2, the dynamics of P-TEGs with strict initial conditions are not pure LDIs. In this subsection, we prove that they can be expressed by means of SLDIs under specific types of schedules, with the immediate consequence that any property of SLDIs also holds for P-TEGs with strict initial conditions. We want to prove that inequalities can be written as SLDIs. For this aim, let us define an auxiliary variable x(0) ∈ R n . Note that the first inequality of (8) is equivalent to where E ij = 0 for all i, j ∈ 1, n . Indeed, the first of the latter inequalities can be rewritten as ∀i ∈ 1, n , max(x 1 (0), . . . , x n (0)) ≤ x i (0) ≤ min(x 1 (0), . . . , x n (0)), which admits as solution all x(0) that satisfy x i (0) = x j (0) for all i, j ∈ 1, n ; therefore, solutions can be parametrized in t 0 ∈ R as x(0) = t 0ẽ . Recalling that P-TEGs (as well as SLDIs) are time-invariant systems (see Subsection 3.1.2), this proves that (8) is equivalent to SLDIs S = ({i, a}, A 0 , A 1 , B 0 , B 1 ) under schedule w = iaaaa . . . = ia K (i.e., mode i followed by a sequence of K modes a), where The following result is an immediate consequence of this fact.

under strict initial conditions determined by matrices ∆ and ∆ is (boundedly) consistent if and only if schedule w = ia +∞ is (boundedly) consistent for the SLDIs S defined as above.
Although Proposition 4 alone does not answer the question about how to verify (bounded) consistency of P-TEGs with strict initial conditions 14 , it suggests that any result regarding SLDIs under schedules of the form w = ia K holds automatically for P-TEGs with strict initial conditions.
In the following, we provide an example for the application of this fact. Let us study the existence of 1-periodic trajectories for P-TEGs with strict initial conditions. From the above discussion, such trajectories correspond to "ultimately" 1-periodic trajectories of the form for SLDIs S = ({i, a}, A 0 , A 1 , B 0 , B 1 ) under schedule w = ia K . We proceed with a strategy similar to the one seen in Subsection 2.4: substituting formula x(k + 2) = λ k x(2) for all k ∈ 1, K − 1 into (6), we get multiplying by λ −k+1 in the third and fourth inequalities, we obtain By defining the extended dater vectorx = [x ⊤ (1) x ⊤ (2)] ⊤ and using Proposition 2, the inequalities can be rewritten in terms ofx as where Clearly, (9) defines a PIC-NCP, whose solution can be found in time O(n 4 ) using Algorithm 1. The conclusion is that periods of consistent 1-periodic trajectories for P-TEGs under strict initial conditions can be obtained in the same strongly polynomial time complexity as for P-TEGs under loose initial conditions.

Analysis of periodic schedules
In this subsection, we analyze bounded consistency and cycle times of SLDIs when schedule w ∈ Σ * ∪ Σ ω is periodic, i.e., when it can be written as w = v K , K ∈ N ∪ {+∞}, and v = v 1 · · · v V ∈ Σ * is a finite subschedule of length V . We define v-periodic trajectories of period λ ∈ R ≥0 for SLDIs under schedule w = v K as those dater trajectories that, for all k ∈ 1, K − 1 , h ∈ 1, V , satisfy x(V k + h) = λx(V (k − 1) + h); Λ v SLDI (S) denotes the set of all periods (or cycle times) λ for which there exists a consistent v-periodic trajectory. Their relationship with 1-periodic trajectories in P-TEGs is illustrated in the following example.
Example 7 Let us analyze the SLDIs S, with Σ = {a, b, c}, and A 0 z , A 1 z , B 0 z , B 1 z defined as in Example 4; now label z ∈ Σ is to be interpreted as a mode. Thus, for each event k, the dynamics of the SLDIs may switch among those specified by the P-TEGs labeled a, b, and c. We consider periodic schedules (ac) K and (ab) K ; observe that for w = v K , with v ∈ {ac, ab} (i.e., v 1 = a and v 2 = c or v 2 = b), the SLDIs following w can be written as: By definingx(k) = [x ⊤ (2(k − 1) + 1), x ⊤ (2(k − 1) + 2)] ⊤ , the above set of inequalities can be rewritten as LDIs: To see the equivalence of (10) and (11), observe that the second block of (11a) reads

From this transformation, we can easily conclude that S is boundedly consistent under v +∞ if and only if the LDIs with characteristic matrices
v are boundedly consistent, and that all consistent v-periodic trajectories of S coincide with consistent 1-periodic trajectories of the LDIs; hence, using Algorithm 1 we obtain Λ ac . It is worth noting that, although P-TEGs labeled a and b are not boundedly consistent, the SLDIs under schedule (ab) +∞ are. Thus, in general it is not possible to infer bounded consistency of SLDIs under a fixed schedule w solely based on the analysis of each mode appearing in w.
By generalizing the procedure shown in Example 7, we can derive the following proposition through some algebraic manipulations (to set up equivalent LDIs) and applying Theorem 3 (for a formal proof, see Appendix A).

Proposition 5 SLDIs S are boundedly consistent under schedule w = v +∞ if and only if they admit a v-periodic trajectory. Moreover, set
Proposition 5 directly provides an algorithm to compute the minimum and maximum cycle times of SLDIs under a fixed periodic schedule. Indeed, these values come from solving the NCP for the parametric precedence graph G(λP v ⊕ λ −1 I v ⊕ C v ). However, this approach results in a slow (although strongly polynomial time) algorithm when the length of subschedule v is large; indeed, its time complexity is O((V n) 4 ) = O(V 4 n 4 ), as the considered precedence graph has V n nodes.
To speed up the computation of Λ v SLDI (S), we may note the following fact: the longer subschedule v is, the larger is the number of −∞'s compared to real numbers in λP v ⊕ λ −1 I v ⊕ C v . In other words, the matrix becomes sparser and sparser with larger values of V . Moreover, the real entries of the matrix have a recognizable pattern. The following theorem, proven in Appendix B, exploits this observation, achieving time complexity O(V n 3 +n 4 ) for computing the set Λ v SLDI (S). The resulting complexity is thus linear in the length of subschedule v.
Theorem 6 Precedence graph G(λPv ⊕ λ −1 Iv ⊕ Cv) does not contain circuits with positive weight if and only if the following conditions hold: The time complexity is analyzed as follows. Item 1 from Theorem 6 requires to check the existence of circuits with positive weight in V precedence The formulas of Theorem 6 generalize those found in [39, Theorem 5.2] for the throughput evaluation in max-plus automata (or heap models). The formula in [39,Theorem 5.2] can indeed be recovered from Theorem 6 considering the special case when P h = C h = E for all h, i.e., when no upper bound constraints or relations between x i (k) and x j (k) for all i, j, k exist; in order to do so, observe that in this particular case Algorithm 1 simplifies significantly, see [27,Remark 3].

Analysis of intermittently periodic schedules
We conclude this section with the analysis of particular trajectories of SLDIs with intermittently periodic schedules, i.e., schedules that can be factorized in the form where • u 0 , . . . , u q , v 1 , . . . , v q are finite subschedules of lengths U 0 , . . . , U q ∈ N 0 and V 1 , . . . , V q ∈ N, respectively, • 2 ≤ m 1 , . . . , m q−1 < +∞, • m q is either an element of N or +∞; in the second case, u q is the empty string (of length U q = 0) and the schedule is called ultimately periodic. We call u 0 , . . . , u q the transient subschedules and v 1 , . . . , v q the periodic subschedules of the schedule w. Observe that the factorization of an intermittently periodic schedule into transient and periodic subschedules may be not unique; for example, schedule w = aababa can be factorized into w = a(ab) 2 a (u 0 = a, v 1 = ab, u 1 = a, m 1 = 2), into w = aa(ba) 2 (u 0 = aa, v 1 = ba, m 1 = 2), or even into w = (a) 2 (ba) 2 (U 0 = U 1 = U 2 = 0, v 1 = a, v 2 = ba, m 1 = m 2 = 2). In the reminder of the paper, to unequivocally indicate the intended schedule factorization into periodic and transient subschedules, we adopt the practice of writing explicitly exponents m h that elevate periodic subschedules v h , and of writing extensively transient subschedules u h (even when u h could be written as a concatenation of a sequence of modes).
The interpretation of intermittently periodic schedules in systems modeled by SLDIs is as follows: after the start-up of the system (u 0 ), a number of operations are executed cyclically (v m1 1 ), after which the system is re-initialized (u 1 ) before starting a new sequence of cyclical tasks (v m2 2 ), and so on; finally, after a finite number (q) of alternations between transient and cyclic working regimes, the system is either shut down (u q ) if m q ∈ N or works in periodic regime forever if m q = +∞.
In the reminder of the paper, we let K h = U 0 + h j=1 (m j V j + U j ) for all h ∈ 0, q . The objective of this section is to show that trajectories with important practical relevance under intermittently periodic schedules can be efficiently analyzed in SLDIs. Namely, we study the existence of intermittently periodic trajectories, that is, trajectories of the dater function {x(k)} k∈ 1,Kq such that are v h -periodic trajectories of period λ h for all h ∈ 1, q . In other words, for all h ∈ 1, q , j ∈ 1, m h − 1 , k ∈ 1, V h , an intermittently periodic trajectory satisfies, for some λ h ∈ R ≥0 , Intermittently periodic trajectories in which m q = +∞ are referred to as ultimately periodic trajectories. Intermittently periodic trajectories generalize v-periodic trajectories, as they are v h -periodic in each sequence of cyclical tasks v m h h . Note that the definition of intermittently periodic trajectory depends on the specific factorization of schedule w into transient (u h ) and periodic (v h ) subschedules. For instance, let w = aababa. A trajectory {x(k)} k∈ 1,6 for schedule w factorized as a(ab) 2 a is intermittently periodic if it satisfies, for some λ 1 ∈ R ≥0 , x(4) = λ 1 x(2), x(5) = λ 1 x(3). Considering the factorization w = aa(ba) 2 , the trajectory is intermittently periodic if, for some λ 1 ∈ R ≥0 , According to factorization w = (a) 2 (ba) 2 , instead, the trajectory is intermittently periodic if, for some λ 1 , λ 2 ∈ R ≥0 , Trajectories of this type appear frequently in applications. Typical examples are batch manufacturing systems, where each batch is processed in a periodic workflow, but switching between different batches leads to pauses and irregular transients (see, e.g., [40,41]). Urban railway systems also operate on a similar principle, with trains arriving at stations in a periodic manner during peak and off-peak hours, although the period may vary based on the time of day. Moreover, note that, as shown in Subsection 4.2, 1-periodic trajectories of P-TEGs with strict initial conditions correspond to a particular case of ultimately periodic trajectories of SLDIs, in which q = 1 and U 0 = 1. for some m 1 ≥ 2. The inequalities corresponding to schedule w are: for all k ∈ 1, m 1 , To analyze the existence of intermittently periodic trajectories, we substitute in (14), for all k ∈ 1, m 1 − 1 , x(k + 2) = λ k 1 x(2), obtaining: for all k ∈ 1, m 1 , It is possible to get rid of the term λ m1−1 1 in the penultimate inequalities by performing a change of variable. Let ξ : 1, 3 → R 2 be defined by By substituting ξ into (15), we obtain Finally, by rewriting the inequalities in terms of the extended vector [ξ ⊤ (1), ξ ⊤ (2), ξ ⊤ (3)] ⊤ and using Proposition 2, we can conclude that the set of λ 1 's for which the SLDIs admit an intermittent periodic trajectory under schedule w coincides with the solution set ΛNCP(λ 1 P ⊕ λ −1 1 I ⊕ C) of the PIC-NCP defined by matrices z for all z ∈ {a, b, c}. Using Algorithm 1 we can compute ΛNCP(λ 1 P ⊕ λ −1 1 I ⊕ C) = [1,1]. In this simple example, the set of admissible periods coincides with Λ c SLDI (S); at a superficial glance one could think that this is a general fact, i.e., that for any intermittently periodic schedule of the form (12), the set of valid periods λ h defined as in (13) is simply Λ v h SLDI (S), for all h ∈ 1, q . This is however not the case, as shown in the following example.
Example 9 (starving philosophers problem, cont.) We analyze the existence of intermittently and ultimately periodic trajectories for Example 6 under schedule where m 1 ∈ N and m 2 = +∞. For the first m 1 dining cycles, schedule w forces the first philosopher to eat twice in a row simultaneously with the third philosopher, who eats once, after which it is the turn of philosophers 2 and 4 to eat; the subsequent dining cycle is identical to the previous m 1 cycles except that the 1 st philosopher eats only once; thereafter the dining order is the one from Exercise 6. Through computations analogous to those seen in the previous example, we can show that all periods λ 1 and λ 2 (corresponding to the two periodic subschedules v 1 = p 1 p 1 p 3 p 2 p 4 and v 2 = p 2 p 4 p 1 p 3 p 3 in w) of consistent intermittently periodic trajectories are those in ΛNCP(A(λ 1 , λ 2 )) = ΛNCP(λ 1 P 1 ⊕ λ −1 Thus, the existence of intermittently periodic trajectories can be checked by solving an MPIC-NCP. For instance, an intermittently periodic trajectory that minimizes the sum of periods λ 1 and λ 2 can be found by solving the following linear programming problem 15 : min The constraints of the above problem can be expressed as 212 inequalities (with each inequality corresponding to an element of A(λ 1 , λ 2 ) different from −∞) in 77 variables. Clearly, the optimization problem can be solved efficiently by linear programming solvers such as interior point methods or the simplex method. However, since the number of optimization variables and constraints increases (linearly) with the length of subschedules u 0 , . . . , uq and v 1 , . . . , vq, this method can be impractical for larger problems. This becomes especially evident when considering hard scheduling problems such as the following one: find the schedule w, from a prescribed set of intermittently periodic schedules, which admits the intermittently periodic trajectory with minimal performance index q i=1 λ i . The intrinsic NP-hardness of the problem requires the use of enumeration methods for its exact solution, and relying on standard solutions approaches for problems like (17) to evaluate the performance index can result in a slow optimization procedure. Therefore, it is natural to ask whether an alternative, less expensive technique to solve problems of the form (17) exploiting the structure of the matrix in (16) exists, similarly to what seen in Theorem 6 for periodic schedules; an affirmative answer to this question will be provided in Theorem 8.
Returning to (17), its solution reveals that the optimal value of λ 1 ⊗ λ 2 is 19, corresponding to λ 1 = 11 and λ 2 = 8. The Gantt chart of Figure 11 shows a trajectory corresponding to a solution of the optimization problem, for m 1 = 2.  Figure 11: Gantt chart representing an intermittently periodic trajectory for the starving philosophers problem.
Before concluding the example note that, for the SLDIs S modeling the starving philosophers problem, Λ p2p4p1p3p3 SLDI (S) = [7.5, 16] (the lower bound is the period of the trajectory shown in the Gantt chart of Figure 10) whereas the value of λ 2 solving (17) is 8 > 7.5. This shows a remarkable fact that is exclusive to systems with upper bound constraints: the cycle times valid in intermittently periodic trajectories may be different from those obtained by considering periodic subschedules independently. The reason for this phenomenon has to be searched in the structure of the precedence graph of matrix A(λ 1 , λ 2 ) (schematically illustrated in Figure 16 on page 51). Here the critical circuit (i.e., the circuit with largest weight) can be generated from arcs connecting portions of the graphs related to different regimes (periodic or transient) of the schedule.
The latter two examples suggested the following proposition, which is proven in Appendix C.

Proposition 7
The set of periods (λ 1 , . . . , λq) ∈ R q ≥0 of intermittently periodic trajectories that are consistent for a given intermittently periodic schedule coincides with ΛNCP( q h=1 (λ h P h ⊕ λ −1 h I h ) ⊕ C), where P h , I h , C are appropriately defined square matrices with (U 0 + q h=1 V h + U h )n rows and columns.
Similarly to Theorem 6, it is possible to reduce the complexity of the problem by leveraging the sparsity and modularity of matrix With the term modularity, we mean the recognizable block structure of the matrix; for an illustrative example see the dashed and dotted blocks in matrix A(λ 1 , λ 2 ) of (16). The main result, which is based on an extensive application of Theorem 6, is stated below. Details of this result and its proof is provided for the particular case of Example 9 in Appendix D.
Theorem 8 The MPIC-NCP of Proposition 7 can be transformed into an equivalent one where the precedence graph to be studied has qn nodes (instead of (U 0 + q h=1 V h + U h )n). The reduction requires O U 0 + q h=1 V h + U h n 3 operations.
Observe that performing the reduction takes a number of operations that is linear in the sum of the lengths of subschedules u 0 , . . . , u q and v 1 , . . . , v q . As no linear programming solver with linear worst-case complexity has ever been found, the advantage of the reduction becomes more prominent with longer subschedules.
Example 10 (starving philosophers problem, cont.) Let A(λ 1 , λ 2 ) be as in (16). Applying Theorem 8 as explained in Appendix D, it can be shown that G (A(λ 1 , λ 2 )) ∈ Γ if and only if G ( A(λ 1 , λ 2 )) ∈ Γ, where A(λ 1 , λ 2 ) = λ 1 P 1 ⊕ λ −1 1 I 1 ⊕ λ 2 P 2 ⊕ λ −1 2 I 2 ⊕ C, and P 1 , I 1 , P 2 , I 2 , C are matrices of dimension 10 × 10 with coefficients in Rmax. In particular, the linear programming problem in (17) has the same optimal value as the following one: min It can be verified that the constraints of the above problem can be expressed as 146 inequalities in 12 variables. Hence, compared with the constraints in (17), we have a 31% reduction in the number of inequalities and 84% reduction in the number of variables.

Practically motivated example
The example we present is a multi-product processing network taken from [26]. Examples of such networks are electroplating lines and cluster tools. Consider a manufacturing system consisting of 5 processing stations S 1 , . . . , S 5 and a robot of capacity one. The system can treat two types of parts, part a, which requires to be processed on S 1 , S 3 , and S 5 (in this order), and part b, which must follow route S 2 , S 1 , S 4 , S 5 . The task of the robot is to transport parts of type a and b from an input storage S 0 to their first processing stations, between the processing stations (in the right order), and finally from the last processing station to an output storage S 6 . The time the robot needs to travel from S i to S j is τ ij when it is not carrying any part, and τ z ij when it is carrying part z ∈ {a, b}. Moreover, the processing time for part z in station S i must be within the interval ι z i = [L z i , R z i ] ⊂ R ≥0 . We consider the following parameters for the processing network:

Cycle time analysis
A classical scheduling problem in this kind of processing networks is to find a periodic robot operation sequence that minimizes the cycle time and avoids time-window constraints violations; in the literature, this is referred to as the single-robot (or hoist) cyclic scheduling problem [26]. Such an optimization problem, which is strongly NP-hard, can be divided into two subproblems: P1 the cycle time minimization, given a fixed sequence of robot operations, P2 the search for the optimal robot operations sequence.
In this section we focus on P1, which is polynomial-time solvable: clearly, an efficient algorithm for P1 can be used as subroutine in search procedures to solve the more general scheduling problem (P1+P2). To find the cycle time for a given robot operation sequence, we will firstly assume, as is standard in robotic cyclic scheduling problems [26,42,43], that the system is already treating parts in a periodic manner from the initial time; the application of Theorem 6 will then provide the possible periods of such treatment plans. In particular, we suppose that initially station S 3 is processing a part of type a, and S 2 , S 4 are processing parts of type b. Of course, this assumption is not met by real systems at start-up time (when all stations are empty), and will thus be relaxed in the next section.
We denote by S i z − → S j robot operation "unload a part of type z from S i , transport it to and load it into S j " and by → S j operation "travel from the current location to S j and wait if necessary". A schedule for this process is a sequence of modes w ∈ {a, b} ω , where mode a represents the sequence of operations Initially, the robot is positioned at S 3 if w 1 = a or at S 4 if w 1 = b. Let us first model the processing network when only part a, respectively, b is considered. In this way, we obtain two P-TEGs, P-TEG a and P-TEG b (shown in Figure 12), each of which represents the behavior of the system when processing only parts of one type. Using Algorithm 1, we can find that the cycle times of the network when processing only parts of type a, b are all values in [73, +∞[, and [72, 192], respectively. Now, from the obtained P-TEGs, we can model the processing network in the case where both parttypes are considered as SLDIs S = ({a, b}, A 0 , A 1 , B 0 , B 1 ). To do so, we must define matrices A 0 z , A 1 z , B 0 z , B 1 z ∈ R n×n max for z ∈ {a, b} appropriately: we start by adding in P-TEG a (respectively, P-TEG b ) the missing transitions from P-TEG b (respectively, P-TEG a ) -the obtained P-TEGs have both n = 12 transitions (in general, n = 2 + 2 × number of processing stations). For each new transition t i of P-TEG z , we define (A 1 z ) ii = (B 1 z ) ii = 0; this is done to store in auxiliary variables x i (k) the last entrance and exit times of parts in stations that are not used in mode z. Moreover, to model the transportation of the robot from S 3 to S 4 (respectively, from S 4 to S 3 ) after each switching of mode from a to b (respectively, from b to a), we set (A 1 a ) 4out,3in = τ 34 (respectively, (A 1 b ) 3out,4in = τ 43 ). The other elements of A 0 z , A 1 z , B 0 z , B 1 z are taken from the characteristic matrices of P-TEG z , for z ∈ {a, b}. This modeling procedure (a) P-TEGa for parts of type a. t0 (b) P-TEG b for parts of type b. Figure 12: P-TEGs modeling the processing network considering only types of one part. A token in a place colored red and black, respectively blue, represents a part being processed in a station and the robot moving with, respectively without carrying a part.
results in the following matrices for mode a: The modeling effort required to define S is repaid by the possibility to use Theorem 6 for computing the minimum and maximum cycle times corresponding to a schedule w = v K , for K ∈ N ∪ {+∞}. For instance, we get Λ ba SLDI (S) = [77,192]. This means that, using schedule (ba) ω , the time between subsequent completions of a product of the same type is at least 72 and at most 192 time units.
To appreciate the advantage of using the algorithm derived from Theorem 6, in Figure 13 we compare the computational time to obtain Λ v SLDI (S) with increasing subschedule length V using Theorem 6 and other methods; specifically, we also consider the algorithm derived from Proposition 5 directly, the algorithm developed in [26], and a linear programming solver. The implementations were made on Matlab R2019a, and for solving the linear programs we used CPLEX's dual simplex method; the tests were executed on a PC with an Intel i7 processor at 2.20Ghz. From the results, we can see that the most time-consuming approach is the one using Proposition 5 directly, while the algorithm from Theorem 6 achieves the fastest computation. The advantage becomes more evident with larger subschedule lengths: for instance, when V = 300, the dual simplex method takes 11.0 · 10 −2 seconds to solve the problem, whereas the algorithm derived from Theorem 6 only 1.25 · 10 −2 seconds. Such computational time reduction may have considerable impact for the solution of cyclic scheduling problems.

Considering start-up and shut-down transients
At the start-up, the system cannot follow the periodic trajectories found in the previous section as all stations are initially empty. Moreover, the periodic workflow must be interrupted for the system shut-down, at the end of which all stations are left empty. In the following, we also suppose that the initial position of the robot coincides with the location of the input storage S 0 .
To represent the complete dynamics of the processing network, from the start-up to the shut-down, we introduce additional modes of operations modeling the initial and final transients. In particular, we add three modes for the start-up: i b1 corresponds to the sequence of operations The subschedule i b1 i b2 i a represents the initial transient of the processing network, consisting of the transportation inside the system of the first two parts of type b and the first part of type a; thus, at the end of the sequence of operations corresponding to i b1 i b2 i a , a part of type b is in station S 2 , another part of the same type is in S 4 , and a part of type a is in S 3 . Similarly, three modes are added for the shut-down: mode f b1 corresponds to Similarly to modes a and b, we can derive matrices A 0 z , A 1 z , B 0 z , and B 1 z for the additional modes i b1 , i b2 , i a , f b1 , f a , f b2 .
An example of complete schedule for the processing network is the intermittently periodic schedule where m ∈ N is the number of repetitions of subschedule ba. To find the admissible periods λ of an intermittently periodic trajectory following schedule w, we solve the PIC-NCP on matrix using Theorem 8; the outcome is that the valid periods are those in interval [77,192]. Finally, a consistent intermittently periodic trajectory can be obtained from any column of (λP ⊕ λ −1 I ⊕ C) * for λ ∈ [77,192] (see Proposition 1). For instance, the trajectory displayed in Figure 14 is derived through some simple manipulations from the first column of matrix A = (77P ⊕ 77 −1 I ⊕ C) * ; in particular, we consider w = i b1 i b2 i a babaf b1 f a f b2 (i.e., m = 2) and impose: In the following, we draw comparisons between SLDIs and other related dynamical systems to provide a broader context for this new class of systems and to outline possible research directions inspired from previous work.
SLDIs have strong connections with several other dynamical systems, with interval weighted automata being the most closely related in terms of modeling expressiveness [44,45]. Interval weighted automata represent the natural extension of max-plus automata to the case of time-window constraints, which force the dater function to satisfy inequalities of the form Expanding the seminal work [39] of Gaubert and Mairesse, in [45] it was shown that safe P-time Petri nets, in which the number of tokens per place cannot exceed 1, are a subclass of interval weighted automata. Comparing (6) with (19), it is not difficult to see that SLDIs are an even larger class of systems compared to interval weighted automata, implying that also SLDIs can represent the dynamics of safe P-time Petri nets. If from (19) we eliminate the upper bound constraints (by defining B 1 w k = T ) and take the least consistent trajectory of the dater function (by substituting the left "≤" sign in (19) by "="), we get the dynamics of a max-plus automaton [46]: This shows that the behavior of any max-plus automaton corresponds to the "fastest" trajectory of specific switched max-plus linear-dual inequalities, in which A 0 w k = E and B 0 w k = B 1 w k = T . From this relationship one can observe that the algorithm derived in this paper for the cycle time computation in periodic schedules (Theorem 6) generalizes [39,Theorem 5.2], which gives a simple formula for the cycle time of safe job shops (without upper bound constraints). The involvement of more complex formulas featured in Theorem 6 has to be attributed to the greater number of circuits in the precedence graph when accounting for upper-bound constraints (see Figure 15; the case considered in [39] corresponds to the same type of graphs but without arcs labeled P i , λP i , and C i ).
In [47], van den Boom and De Schutter extended the capabilities of maxplus automata by introducing an input signal u(k) ∈ R m and an input-state matrix D w k ∈ R n×m max to (20); the resulting dynamical system is referred to as a switching max-plus linear system: Now the system is not anymore forced to evolve according to the fastest possible trajectory represented by (20), but it is possible to construct a controller that, by selecting the appropriate input signal, can delay the occurrence of events with the aim of regulating the behavior of the system. This feature makes switching max-plus linear systems, in a certain sense, more similar to SLDIs. We explain this in more detail in the following. Take a switching max-plus linear system and assume that the input signal is able to delay "directly" the occurrence of each event. Mathematically, this corresponds to having D w k = E ⊗ ∈ R n×n max for all k. Then, recalling that ⊕ in (21) is the elementwise maximization, by taking u(k + 1) sufficiently large (in particular, A 1 w k ⊗ x(k) ≤ u(k + 1)), x(k + 1) and u(k + 1) coincide, and (21) simplifies to which is equivalent to A 1 w k ⊗ x(k) ≤ x(k + 1). In this way, we have recovered a portion of the inequalities that appear in SLDIs. Thus, observe that: • on the one hand, unlike SLDIs, switching max-plus linear systems are not able to represent upper-bound constraints, • on the other hand, only systems where all events are directly controllable (i.e., delayable) can be modeled by SLDIs, whereas this does not need to be the case for switching max-plus linear systems. It follows that the modeling expressiveness of switching max-plus linear systems and SLDIs is not comparable. The generalization of the results of the present paper to switching dynamical systems driven by time-window constraints and where the assumption of direct controllability of all events is relaxed is left as future work.
Regarding control approaches, it would be valuable to extend to SLDIs techniques already investigated in switching max-plus linear systems and maxplus automata: we mention model predictive control [47], just-in-time control based on residuation theory [48], geometric control [49], and supervisory control [50].
To conclude, we emphasize that interesting connections can be drawn between the concepts of stability in discrete-time switched linear systems 16 (see, e.g., [51]) and bounded consistency in SLDIs. Considering the case of "unswitched" systems (i.e., discrete-time linear systems and LDIs), the two properties are both verifiable in strongly-polynomial time (using Jury's test for the former property [52], using Theorem 3 for the latter) and are related to the spectral radius in the respective algebras. In their switching counterparts, stability and bounded consistency are -unsurprisingly -interconnected with the notion of joint spectral radius in the standard and max-plus algebras. Further investigation of this link presents an exciting avenue for future research. 16 Recall that discrete-time switched linear systems follow x(k + 1) = A 1 w k · x(k), where · indicates the standard matrix-vector multiplication.

A Proof of Proposition 5
Let us consider SLDIs S under schedule By definingx and using Proposition 2, it is possible to rewrite the SLDIs above as the following LDIs: for all k ∈ N, Hence, since SLDIs under periodic schedule v +∞ are equivalent to the above LDIs, according to Theorem 3 they are boundedly consistent if and only if they admit, for some λ ∈ R, a trajectory of the form ∀k ∈ N,x(k) = λ k−1x (1); observe that such a trajectory is v-periodic as, in terms of dater x, it satisfies the following property: To find the set Λ v SLDI (S) of admissible periods for the SLDIs under schedule v +∞ , we proceed as in Section 2.4, obtaining: Therefore, the set Λ v SLDI (S) coincides with all the λ's such that the precedence

B Proof of Theorem 6
In this section, we prove Theorem 6. All symbols used in the following propositions are as in Proposition 5 and Theorem 6. To prove the theorem, we first need a number of technical propositions and lemmas, which provide closed formulas for the computation of elements of the Kleene star of particular sparse matrices. We start from recalling some properties of the Kleene star operation in the max-plus algebra.

Proposition 9 ([22])
Given two matrices a, b ∈ R n×n max , the following properties hold: Some of the following intermediate results have already been concisely proven in [53]; here we provide a more detailed proof for the first of them.
The lemma is proven by applying the formulas from the latter proposition to Algorithm 2 from [24], which can be used for the computation of the Kleene star.
where each M ij is a matrix partitioning M * such that M 11 ∈ R n1×n1 max ; applying Algorithm 2 of [24], we easily get = a * bd * (ca * bd * ) * ca * ⊕ a * = a * a * bd * d * (ca * a * bd * d * ) * ca * a * ⊕ a * a * = a * a * bd * (d * ca * a * bd * ) * d * ca * a * ⊕ a * a * where a + = aa * and, consequently, a * = a + ⊕ E ⊗ . The formula for matrix M 12 can be obtained in a similar way starting from which comes from Algorithm 2 of [24]. Formulas for M 21 and M 22 are directly derived from reasoning about the symmetry of the Kleene star operation.
Let us now divide matrix λP v ⊕ λ −1 I v ⊕ C v into four blocks, such that the top-left block is C 1 ∈ R n×n max . We will first focus on formulas to compute specific elements of the Kleene star of the bottom-right block of λP v ⊕ λ −1 I v ⊕ C v , which we will denote by M . We use the following notation: where each M ij is an n × n matrix. Note that matrix M is a block tridiagonal matrix, i.e., a matrix whose blocks are all E except for those in the main diagonal and in the first diagonals above and below the main diagonal. We refer to [54,55] for some work related to tridiagonal matrices in the tropical (i.e., either max-plus or min-plus) algebra.
Proof The proof is done by induction on V . For V = 3, the formula can be obtained directly from Lemma 10: M 32 = (C * 3 I 2 C * 2 P 2 C * 3 ) * C * 3 I 2 C * 2 = (C I 2 ) * I 2 . Suppose now that the formula holds for V −1 with V ≥ 4; we prove that it holds also for V . By partitioning matrix M in four blocks, such that the bottom-right block coincides with C V , from Lemma 10 we get that M V 2 is equal to the leftmost block of matrix  Expanding the expression using the induction hypothesis, Lemma 12,(23), and (24), we can get the desired formula with computations analogous to those in the proof of Lemma 13.

Remark 3
With similar techniques, one can prove that, for all V ≥ 4 and for all i ∈ 3, V − 1 , In the next proposition, we focus on the top-left n × n block of matrix (λP v ⊕ λ −1 I v ⊕ C v ) * ; we denote this block by N ∈ R n×n .
Proof Partitioning matrix λPv ⊕ λ −1 Iv ⊕ Cv into four blocks, such that the top-left block is C 1 , Lemma 10 gives us where the latter step can be obtained from (23), and (24). Finally, the commutativity of ⊕ concludes the proof.
After the technical propositions and lemmas above, we are finally ready to introduce the graphical part of the proof. First of all, it is useful to visualize the structure of G(λP v ⊕λ −1 I v ⊕C v ); the lumped-node representation of Figure 15  illustrates it in the case V = 5. In this simplified representation, k indicates the set of nodes (k − 1)n + 1, kn for any k ∈ 1, 5 , and the matrix, say Y , associated to an arc from k 1 to k 2 indicates that an arc in G(λP v ⊕λ −1 I v ⊕C v ) from node (k 1 − 1)n + j to node (k 2 − 1)n + i exists iff Y ij = −∞, and that its weight is equal to Y ij , for all i, j ∈ 1, n .
Recall that, for a square matrix M such that M * = M 11 M 12 M 21 M 22 , where M ii has dimension n i × n i , G(M ) does not contain circuits with positive weight visiting any node j ∈ 1, n 1 (respectively, j ∈ n 1 + 1, n 2 ) iff G(M 11 ) ∈ Γ (respectively, G(M 22 ) ∈ Γ). We partition the set of circuits of G(λP v ⊕λ −1 I v ⊕ C v ) in V subsets as follows: (1) circuits visiting at least one node in 1, (2) circuits that do not visit any node in 1, but that visit at least one node in 2, (3) circuits that do not visit any node in 1 or 2, but that visit at least one node in 3, . . . , (V ) circuits that do not visit any node in 1, 2, . . . , V − 1, but that visit at least one node in V . From Proposition 15, the maximal weights of circuits in the first subset correspond to elements in the diagonal of N = C * 1 (λM P ⊕ λ −1 M I ⊕ M C ) * C * 1 ; thus these are non-positive iff G(C 1 ) ∈ Γ and λ ∈ Λ NCP (λM P ⊕ λ −1 M I ⊕ M C ). Let k ∈ 2, V − 1 . From Lemma 11, the maximal weights of circuits in the k th subset correspond to elements in the diagonal of C * k (C P k ) * C * k ; thus, these are non-positive iff G(C k ) ∈ Γ and G(C P k ) ∈ Γ. As for circuits in the V th subset, their maximal weights come from the diagonal elements of C * V , which are finite iff G(C V ) ∈ Γ.

D Proof of Theorem 8
For notational simplicity, we show how to prove the theorem for the matrix in (16); the extension to the general case is simple but cumbersome.
First of all, it is useful to draw the parametric precedence graph G(A(λ 1 , λ 2 )) using the lumped-node representation (which was presented in Appendix B); the graph is shown in Figure 16. We recall that the objective is to describe in a more compact way the set Λ NCP (A(λ 1 , λ 2 )) of all parameters (λ 1 , λ 2 ) for which G(A(λ 1 , λ 2 )) does not contain circuits with positive weight. Before illustrating the proof, it is convenient to relabel the nodes of G(A(λ 1 , λ 2 )) as in Figure 17; the graph in Figure 17 therefore corresponds to the precedence graph G(A(λ 1 , λ 2 )), where A is obtained permuting rows and columns of A(λ 1 , λ 2 ) as follows:  Figure 17: Lumped-node representation of G(A(λ 1 , λ 2 )) from (38).
To check the second condition, we must compute a * bd * ca * . The main challenge here is to find a closed formula for Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2002/1 "Science of Intelligence" -project number 390523135.