Synchronization of a class of cyclic discrete-event systems describing legged locomotion

It has been shown that max-plus linear systems are well suited for applications in synchronization and scheduling, such as the generation of train timetables, manufacturing, or traffic. In this paper we show that the same is true for multi-legged locomotion. In this framework, the max-plus eigenvalue of the system matrix represents the total cycle time, whereas the max-plus eigenvector dictates the steady-state behavior. Uniqueness of the eigenstructure also indicates uniqueness of the resulting behavior. For the particular case of legged locomotion, the movement of each leg is abstracted to two-state circuits: swing and stance (leg in flight and on the ground, respectively). The generation of a gait (a manner of walking) for a multiple legged robot is then achieved by synchronizing the multiple discrete-event cycles via the max-plus framework. By construction, different gaits and gait parameters can be safely interleaved by using different system matrices. In this paper we address both the transient and steady-state behavior for a class of gaits by presenting closed-form expressions for the max-plus eigenvalue and max-plus eigenvector of the system matrix and the coupling time. The significance of this result is in showing guaranteed robustness to perturbations and gait switching, and also a systematic methodology for synthesizing controllers that allow for legged robots to change rhythms fast.


I. INTRODUCTION
Synchronization of cyclic processes is important in many fields, including manufacturing [51], transportation [31], genomics [48], and neuroscience [34,50], etc (see references within [13]). In this paper we focus on a class of multiple concurrent two-state cyclic systems with a direct application to legged locomotion. Our motivation is the requirement of legged mobile robots to cope with unstructured terrains in a flexible way by smoothly and effectively switching between different gaits.
Legged systems are traditionally modeled using cross-products of circles in the phase space of the set of continuous time gaits. Holmes et al. [34] give an extensive review of dynamic legged locomotion.
The central pattern generator [35] approach to design motion controllers lies in assembling sets of error functions to be minimized that cross-relate the phases of multiple legs, resulting in attractive limit cycles for the desired gait. Switching gaits online is typically not addressed since in the central pattern generator framework switching must be modeled as a hybrid system. Additionally, implementing "hard constraints" on the configuration space [29] can be quite complex mainly due to the combinatorial nature of the gait space, and often it comes at the cost of dramatically increasing the complexity of the controller.
As an alternative to the common continuous time modeling approach, we introduce an abstraction to represent the combinatorial nature of the gait space for multi-legged robots into ordered sets of leg index numbers. This abstraction combined with max-plus linear equations allows for systematic synthesis and implementation of motion controllers for multi-legged robots where gait switching is natural and the translation to continuous-time motion controllers is straightforward [38]. The methodology presented is particularly relevant for robots with four, six, or higher numbers of legs where the combinatorial nature of the gait space starts to play an important role. For a large number of legs it is not obvious in which order each leg should be in swing or in stance. Most legged animals, in particular large mammals, are known to walk and run with various gaits on a daily basis, depending on the terrain or how fast they need to move. The discrete-event framework presented in this paper enables the same behavior for multi-legged robots. Mathematical properties are derived for this framework, giving extra insight into the resulting robot motion.
The main mathematical representation employed in this paper are switching max-plus linear equations.
By defining the state variables to represent the time at which events occur, systems of linear equations in the max-plus algebra [1,11,32] can model a class of timed discrete-event systems. Specifically, max-plus liner systems are equivalent to (timed) Petri nets [44] where all places have a single incoming and a single outgoing arc. Max-plus linear systems inherit a large set of analysis and control synthesis tools thanks to many parallels between the max-plus-linear systems theory and the traditional linear systems theory. Discrete-event systems that enforce synchronization can be modeled in this framework. Max-plus algebra has been successfully applied to railroads [3,31], queuing systems [30], resource allocation [19], and recently to image processing [2] and legged locomotion [38,39].
The contributions of this paper are the following: we present a class of max-plus linear systems that realize the synchronization of the legs of a robot. Next, we derive closed-form expressions for the maxplus eigenvalue and eigenvector of the system matrix, and show that the max-plus eigen-parameters are max-plus unique, implying a unique steady-state behavior. This result is then used to compute the coupling Existence of eigenvalue and eigenvector (algebraic computation, Lemma 11) Coupling time (algebraic computation, Lemma 16) Transient Irreducibility (matrix structure, Lemma 12) Unique eigenvalue (Corollary 13) Single strongly connected subgraph (graph approach, Unique eigenvector (Theorem 15) Existence Uniqueness Steady-state time, which characterizes the transient behavior. The importance of having closed-form expressions and uniqueness of the max-plus eigen-structure is that, not only can one compute these parameters very fast without recurring to simulations or numerical algorithms (e.g. Karp's algorithm [1]), but one has also guarantees of uniqueness: the motion of the robot will always converge in a finite number of steps to the same prescribed behavior, regardless of gait changes or disturbances. This reassurance is fundamental when designing gait controllers for robotics. Additionally we present a low least number of steps needed to reach steady-state motion after changing gaits. This paper is focused on the general mathematical properties of a class of discrete-event systems that describe legged locomotion, and not on the actual implementation of gait controllers for robots, as presented previously in [38].
In Section II we revisit the fundamentals of legged locomotion with special emphasis on gait generation and show that max-plus algebra can be used in the modeling of the synchronization of multiple legs. In Section III we briefly review relevant concepts from the theory of max-plus algebra and in Section IV we present a class of parameterizations for the gait space. Given such class, we derive a number of properties, such as the max-plus eigen-structure of the system matrices (Section IV-C), their graph representation (Section IV-D), and the coupling time (Section IV-E). Figure 1 illustrates the structure of the contributions presented in this paper.
Modeling of legged locomotion. a) The configuration space of two oscillators is a torus. Synchronization is achieved by constructing attractive limit cycles, represented by the curve on the torus. b) Discrete-event representations of multiple circuits can be modeled as concurrent cyclic Petri nets. Synchronization is achieved by adding extra places, represented by s2 and s2.

II. MODELING LEGGED LOCOMOTION
In literature [14,28,34,45,46], the study of legged locomotion is approached from two main directions: the signal generation side, where emphasis is placed on the classes of signals that result in periodic locomotion behavior independent of the physical platform; and the mechanics side, where the (hybrid) Newtonian mechanics models are analyzed independently of the driving control signal. We focus on the first approach, by restating the traditional view of periodic gaits for legged systems being defined in the n-torus: Cartesian products of circles each representing an abstract phase that parameterizes the position of each leg in the Euclidean space. Such an abstraction serves as a platform for the models of "networks of phase oscillators" and central pattern generators, introduced in the earlier works of Grillner [27] and Cohen et al. [5]. These are now accepted by both biology and robotics communities as standard modeling tools [34].
In this paper we abstract beyond the notion of a continuous phase to consider concurrent circuits of discrete events. Taking a Petri net modeling approach, during ground locomotion the places represent leg stance (when the foot is touching the ground and supporting the body) and leg swing (when the foot and all parts of the leg are in the air). The transitions represent leg touchdown and lift off. This labeling is most convenient for ground locomotion, but one should be aware that the framework presented in this paper is valid for other types of locomotion where the phases of the various limbs need to be synchronized, such as swimming or flight. Figure 2 illustrates the conceptual difference between our approach ( Figure   2.b) where each phase is represented by a circuit in the discrete-event systems domain and the traditional continuous phase central pattern generator (Figure 2.a). Here, phase synchronization is enforced on the torus by implementing controllers that achieve stable limit cycles [37]. In this paper we translate timed event graphs 1 into the equivalent representation as max-plus linear systems, and achieve synchronization 1 We restrict ourselves to a class of timed Petri nets called timed event graphs such that a one-to-one translation to max-plus linear systems is possible (see [32], chapter 7). In timed event graphs each place can have one single incoming arc and one single outgoing arc. . 3. Illustration of a walking pattern, photos by Muybridge [40]. The solid bars, following Hildebrand's diagram notation [33], indicate that the leg is in stance (foot touching the ground) for τg time units, and white space that the leg is in swing (foot in flight) for τ f time units.
The time length when both feet are touching the ground is called the double-stance time τ∆. We use the notation t1(k − 1) to represent the time instant when leg 1 (left leg) touches the ground, and l1(k) when it initiates a swing. The parameter k is the "step counter". by designing the system matrices appropriately. Consider the leg synchronization of a biped robot. In this case, only two legs need to be synchronized during locomotion. For a typical walking motion (no aerial phase) the left leg should only lift off the ground after the right leg has touched down, to make sure the robot does not fall due to lack of support. This simple synchronization requirement can be captured by introducing state variables for the transition events defined as follows: let l i (k) be the time instant leg i lifts off the ground and t i (k) be the time instant it touches the ground, both for k-th iteration, where k is considered to be a global event or "step" counter. Enforcing that the time instant when the leg touches the ground must equal the time instant it lifted off the ground for the last time plus the time it is in swing A similar relation can be derived for the lift off time: where τ g is the stance time and t i (k − 1) refers to the previous iteration such that equations (1) and (2) can be used iteratively. For this system we have that τ f > 0 and also τ g > 0. Synchronization of the cycles of two legs can be achieved by introducing a double stance time parameter, denoted τ ∆ , representing that after each leg touchdown both legs must stay in stance for at least τ ∆ time units (see Figure 3). This is captured by the following equations: Equation (3) enforces simultaneously that leg 1 stays at least τ g time units in stance and will only lift off at least τ ∆ time units after leg 2 has touched down. When both conditions are satisfied, lift off takes place.
Equation (4) is analogous for leg 2. Note that the parameters τ f , τ g , and τ ∆ are in fact the minimal swing, stance, and double stance times, respectively, as opposed to the exact times. Note additionally that these equations are non-linear but they rely only on the "max" and "plus" operations. This motivates the use of the theory of the max-plus algebra to find parsimonious discrete-event models for legged locomotion.
Although the previous set of equations were designed for a walking behavior with no aerial phases, in practice they can also be used for running by choosing the double stance time to be negative, i.e. by enforcing that one leg lifts off τ ∆ time units before the other leg touches the ground. However, the results presented in this paper focus on the case when τ ∆ 0 such that stability (in the sense of ensuring a desired minimum number of legs on the ground simultaneously) is guaranteed.

III. MAX-PLUS ALGEBRA
In the early sixties the fact that certain classes of discrete-event systems can be described by models using the operations max and + has been discovered independently by a number of researchers, among whom Cuninghame-Green [9,10] and Giffler [21][22][23]. These discrete-event systems are called maxplus-linear systems since the model that describes their behavior becomes "linear" when formulated it in the max-plus algebra [1,11,32], which has maximization and addition as its basic operations. More specifically, discrete-event systems in which only synchronization and no concurrency or choice occur can be modeled using the operations maximization (corresponding to synchronization: a new operation starts as soon as all preceding operations have been finished) and addition (corresponding to the duration of activities: the finishing time of an operation equals the starting time plus the duration). Some examples of max-plus linear discrete-event systems are production systems, railroad networks, urban traffic networks, queuing systems, and array processors [1,11,32].
An account of the pioneering work of Cuninghame-Green on max-plus system theory has been given in [11]. Related work has been done by Gondran and Minoux [24][25][26]. In the eighties the topic attracted new interest due to the research of Cohen, Dubois, Moller, Quadrat, Viot [6][7][8], Olsder [41][42][43], and Gaubert [15][16][17], which resulted in the publication of [1]. Since then, several other researchers have entered the field. For an historical overview we refer the interested reader to [12,20,32]. In this section we give an introduction to the max-plus algebra. This section is based on [1,11], where a complete overview of the max-plus algebra can be found. The basic operations of the max-plus algebra are maximization and addition, which will be represented by ⊕ and ⊗ respectively x ⊕ y = max(x, y) and x ⊗ y = x + y, for x, y ∈ R max := R ∪ {−∞}. The zero element for ⊕ in R max is ε := −∞ and the unit element for ⊗ is e := 0. The structure (R max , ⊕, ⊗) is called the max-plus algebra [1,11]. The operations ⊕ and ⊗ are called the max-plus-algebraic addition and max-plus-algebraic multiplication respectively since many properties and concepts from linear algebra can be translated to the max-plus algebra by replacing + by ⊕ and × by ⊗.The max-plus algebra is a typical example of a class of algebraic structures called commutative dioids. In a dioid the additive operation ⊕ is associative, commutative, and idempotent, and it has a zero element; the multiplicative operation ⊗ is commutative and associative and it has an identity element; the additive zero element is absorbing for ⊗; and ⊗ is left and right distributive w.r.t. ⊕. Let r ∈ R. The rth max-plus power of x ∈ R is denoted by x ⊗r and corresponds to rx in conventional There is no inverse element for ε since ε is absorbing for ⊗. If r > 0 then ε ⊗r = ε. If r < 0 then ε ⊗r is not defined. In this paper we have ε ⊗0 = 0 by definition. The rules for the order of evaluation of the max-plus operators are similar to those of conventional algebra. So max-plus power has the highest priority, and max-plus multiplication ⊗ has a higher priority than max-plus addition ⊕. Throughout this paper the i, j element The matrix E n is the n by n max-plus identity matrix: We also define the m by n max-plus "one" matrix 1 m×n such that [1] ij = e = 0 for all i, j. If the dimensions of E, E, 1 are omitted in this paper, they should be clear from the context. The basic max-plus-algebraic operations are extended to matrices as follows.
for all i, j. Note the analogy with the definitions of matrix sum and product in conventional linear algebra.
The max-plus product of the scalar α ∈ R max and the matrix A ∈ R m×n max is defined by [α⊗A] ij = α⊗[A] ij for all i, j. The max-plus matrix power of A ∈ R n×n max is defined as follows: A ⊗0 = E n and A ⊗p = A⊗A ⊗p−1 for p 1.
Theorem 1 (see [1], Th 3.17). Consider the following system of linear equations in the max-plus algebra: If A * exists then solves the system of max-plus linear equations (7) Definition 2. The matrix A ∈ R n×n max is called nilpotent if there exists a finite positive integer p 0 such that for all integers p p 0 we have A ⊗p = E.

It is easy to verify that if
defined as an ordered pair (V,A), where V is a set of vertices and A is a set of ordered pairs of vertices.
The elements of A are called arcs. A loop is an arc of the form . When the initial and the final vertex of a path coincide, we have a circuit. An elementary circuit is a circuit in which no vertex appears more than once, except for the initial vertex, which appears exactly twice. A directed . . , n} and if we associate a real number [A] ij with each arc (j, i) ∈ A, then we say that G is a weighted directed graph. We call [A] ij the weight of the arc (j, i). Note that the first subscript of [A] ij corresponds to the final (and not the initial) vertex of the arc (j, i).

Definition 3 (Precedence graph).
Consider A ∈ R n×n max . The precedence graph of A, denoted by G(A), is a weighted directed graph with vertices 1, 2, . . . , n and an arc (j, i) with weight Let A ∈ R n×n max and consider G(A). The weight |p| w of a path p : i 1 → i 2 → . . . → i l+1 is defined as the sum of the weights of the arcs that compose the path: The average weight of a circuit is defined as the weight of the circuit divided by the length of the circuit: |p| w /|p| l .

Definition 4 (Irreducibility). A matrix A ∈ R n×n
max is called irreducible if its precedence graph is strongly connected.
Definition 5 (Max-plus eigenvalue and eigenvector). Let A ∈ R n×n max . If there exist a number λ ∈ R max and a vector v ∈ R n max with v = E n×1 such that A ⊗ v = λ ⊗ v, then we say that λ is a max-plus eigenvalue of A and that v is a corresponding max-plus eigenvector of A.
It can be shown that every square matrix with entries in R max has at least one max-plus eigenvalue (see e.g. [1]). However, in contrast to linear algebra, the number of max-plus eigenvalues of an n by n matrix is in general less than n. If a matrix is irreducible, it has only one max-plus eigenvalue (see e.g. [6]).
Moreover, if v is a max-plus eigenvector of A, then α ⊗ v with α ∈ R is also a max-plus eigenvector of A.
The max-plus eigenvalue has the following graph-theoretic interpretation. Consider A ∈ R n×n max . If λ max is the maximal average weight over all elementary circuits of G(A), then λ max is a max-plus eigenvalue of A. For formulas and algorithms to determine max-plus eigenvalues and eigenvectors the interested reader is referred to [1,4,6,36] and the references cited therein. Every circuit of G(A) with an average weight that is equal to λ max is called a critical circuit. The critical graph G c (A) of the matrix A is the set of all critical circuits. Let N be the set of positive, non-zero integers.
Theorem 6. Let A be an irreducible matrix. Then there exists c ∈ N (the cyclicity of A), λ ∈ R (the unique max-plus eigenvalue of A), and k 0 ∈ N (the coupling time of A) such that Proof: See e.g. [1,6,18].

IV. LEGGED LOCOMOTION VIA MAX-PLUS MODELING
In Section II equations (1), (3), and (4) describe the synchronization constraints between two legs. We can generalize these equations by defining the following vectors for an n-legged robot: Equation (1) is then written as: If one assumes that the synchronization is always enforced on the lift off time of a leg, equations (3) and (4) are written jointly as: where the matrices P and Q encode the synchronization between lift off of a leg related to a touchdown of the current event (as in equation (4)) and a touchdown of the previous event (as in equation (3)), respectively. The rationale behind this particular model is to prevent that a legged platform has too many legs in swing while walking 2 risking falling down. Synchronization constraints are always imposed on legs that are in stance and are about to enter swing: some legs should only swing if others are in stance (equation (14)). Once in swing, legs are never constrained to go into stance (equation (13)).
Equations (13) and (14) are written in state-space form as: Define the matrices Consider the full state x defined as (15) can then be written in simplified notation: Note that additional max-plus identity matrices E are introduced in the diagonal of matrix A 1 . This results in the extra trivial constraints t i (k + 1) t i (k) and l i (k + 1) l i (k), also resulting in the final system matrix (defined in page 13, equation (35)) being irreducible. This is observed later on in Lemma 12.

A. A gait parameterization
Consider a general legged robot where a two-event circuit is associated to each leg. We present a parsimonious representation of a walking gait of a robot by grouping sets of legs and specifying in what order they are allowed to cycle.
i.e., the sets p form a partition of {1, . . . , n}. A gait G is defined as an ordering relation of groups of legs: The gait space is the set of all gaits that satisfy the previous definitions.
By considering that each p contains the indices of a set of legs that are synchronized in phase, the previous ordering relation is interpreted in the following manner: the set of legs indexed by i swings synchronously. As soon as all legs in i touchdown then all legs in i+1 initiate their swing motion. The same is true for m and 1 , closing the cycle. For example, a trotting gait, where diagonal pairs of legs move synchronously, for a quadruped robot as illustrated in Figure 4, is represented by: sideways [49]. With the previous notation we can now derive the matrices P and Q in equation (15): For the trotting gait G trot we obtain: Define the function that transforms a gait into a vector of integers: For a gait G, define the similarity matrixC ∈ R n×n max as: The similarity matrixC is such thatC ⊗C T =C T ⊗C = E. The similarity matrix associated with the trotting gait G trot is:C here written in both max-plus and traditional algebra notation for legibility purposes. The similarity matrix C has the property of "normalizing" the P and Q matrices to a max-plus algebraic lower triangular form P and a max-plus algebraic upper triangular formQ respectively: Taking the previous example of the trotting gait, the normalized matrices take the form which are generated by the normal gait {1, 2} ≺ {3, 4}. Let # i represent the number of elements of the set i . For a general normal gaitḠ with1 i,j = 1 # i ×# j the structure of the matricesP andQ is: From equation (31) it is clear that the matrixP is always max-plus nilpotent, since the upper triangle is max-plus zero. For max-plus powers ofP we obtain: Max-plus nilpotency is invariant to max-plus similarity transformations (e.g. as defined in equations (27), (28)).
Given an arbitrary gait G with associated matrices P , Q, A 0 , and A 1 one can find the normal matrix P which is max-plus nilpotent. From Lemma 9 then P is also max-plus nilpotent. In the beginning of Section IV we have presented the synchronization equations (17) implicitly. However, if A * 0 exists then using equations (7) and (9), system (17) can be transformed into an explicit set of equations.
Lemma 10. A sufficient condition for A * 0 to exist is that the matrix P is nilpotent in the max-plus sense.
Proof: By direct computation, the repetitive products of A 0 can be found to be If P is max-plus nilpotent, then there exists a finite positive integer p 0 such that ∀p p 0 : = E, and therefore the max-plus sum for the computation of A * 0 is finite: Since P is always max-plus nilpotent for gaits generated by expressions (21) and (22), we conclude that A * 0 is well defined. Let A, that we call system matrix, be defined by: Equation (17) can be rewritten as: For an arbitrary gait the internal structure of A can be quite complex. However, the gait G associated to A can be transformed into a normal gait via a similarity transformation. Let The similarity matrix C transforms the system matrix A of an arbitrary gait G into the system matrix A of a normal gaitḠ via the similarity transformationĀ = C ⊗ A ⊗ C T . This can be shown by direct computation: Transforming an arbitrary gait into a normal gait is very useful since, by effectively switching rows and columns in A, one obtains a very structured matrixĀ where analysis is simple. The interpretation of the similarity matrixC is that legs can be renumbered, simplifying algebraic manipulation. Besides max-plus nilpotency, other properties are invariant to similarity transformations: irreducibility is preserved since the graphs of A andĀ are equivalent up to a label renaming. Max-plus eigenvalues and eigenvectors are related by: B. Structure of the system matrixĀ Let: The structure ofĀ can be obtained via a laborious but straightforward set of algebraic manipulations. For an arbitrary gait G we compute the normal gaitḠ via the similarity transformation with the matrix C. By observing the structures ofĀ 0 andĀ 1 (derived fromP andQ) a closed-form solution can be obtained where W = (τ f ⊗P ) * , illustrated in equation (48) on page 16. The matrixW is defined in equation (47) again on page 16. Note that τ f ⊗W ⊕ E = W and W W . An expression forĀ is then obtained: as illustrated by equation (49). One can show that: Since µ ⊗ W W for any µ > 0, and W W , expression (42) simplifies to: Let Equations (50)-(53) illustrate the resulting structure ofĀ written in the system formx(k) =Ā ⊗x(k − 1), (53), the ⊗ operator is omitted in unambiguous locations due to space limitations.

C. Max-plus eigenstructure of the system matrix
Given a parameterization of a gait, it is fundamental to understand whether the system x(k) = A ⊗ x(k − 1) reaches a unique steady state behavior. In robotics this is the equivalent of asking "does the robot walk/run as specified? Is it robust to disturbances?". These questions are answered by analyzing the max-plus eigenstructure of the system matrix: a unique eigenvalue means that the legs have a unique . . .
cycle time, and a unique (up to scaling) eigenvector means that the legs always reach the same motion pattern, independently of the initial condition or disturbances. The results obtained below use various analysis techniques available for max-plus linear systems. This is necessary due to the intrinsic time structure associated with the problem. Petri net tools (e.g. incidence matrices) can be used to understand structural properties of the system, such as irreducibility, but temporal properties are better analyzed using max-plus linear tools. In Section IV-D we show that for a fixed structure (i.e. a single Petri net) unique or non-unique eigenvectors are found by changing the holding time parameters. This result could not be captured by the Petri net structure alone. The analysis steps presented from here on are summarized in Figure 1.
Consider the following assumption (which is always satisfied in practice since the leg swing and stance times are always positive numbers): Lemma 11. If assumption A1 is satisfied then is a max-plus eigenvalue of the system matrix A (andĀ) defined by equations (36) (and (50)-(53)), and v ∈ R 2n max defined by ∀j ∈ {1, . . . , m}, ∀q ∈ j : is a max-plus eigenvector of A.
Recall equations (7) and (9) with new variables z and B such that Given the previous result, it is sufficient to show that if λ and v are a max-plus eigenvalue and eigenvector of A respectively, then replacing the state variable x(k − 1) by v and x(k) by λ ⊗ v in equation (15) holds true: The previous expression is equivalent to the following two equations: Since λ > 0 (by assumption A1), equation (59) is always verified. Thus we focus on equation (60), which can be simplified due to τ f ⊗ τ g > 0: Let τ ∆ ⊗ P 0 = P and τ ∆ ⊗ Q 0 = Q, i.e., all entries of matrices P 0 and Q 0 are either e or ε to obtain (recall that τ δ = τ f ⊗ τ ∆ and τ γ = τ f ⊗ τ g ): We now consider two cases: i) First we analyze the row indices of equation (62) that are elements of the sets 2 , . . . , m . For each j ∈ {1, . . . , m − 1} and for each row p ∈ j+1 we obtain (notice that according to (22) all the elements of [Q 0 ] p,· are ε since p / ∈ 1 , and that [v] p = τ ⊗j δ for p ∈ j+1 ): The last term always holds true since λ τ γ . Thus for rows p ∈ 2 , . . . , m equation (62) holds true.
ii) We now look at all the remaining rows p such that p ∈ 1 (noticing now that according to (21) all the elements of [P 0 ] p,· are ε and that [v] p = e since p ∈ 1 ): Combining i) and ii) we conclude that equation (62) holds true.
Consider again the trotting gait for a quadruped G trot defined in (20). For this gait m = 2, resulting in: conclude that A is also irreducible.
The max-plus eigenvector v defined by (55)-(56) is not necessarily unique, given assumption A1 alone.
Since, to the authors' best knowledge, there exists no algebraic method to prove max-plus eigenvector uniqueness in general, we take advantage of the precedence graph ofĀ to further investigate this property.
If the critical graph of a irreducible max-plus system matrix has a single strongly connected subgraph, then its max-plus eigenvector is unique up to a max-plus scaling factor (see [1], Theorem 3.101). We proceed by computing the critical graph(s) ofĀ.

D. The precedence graph ofĀ
Given expression (53) it is possible to construct the precedence graph ofĀ. Since this graph can be quite large for a generalĀ, we find it more efficient to first group "similar" nodes into a single node, i.e. apply a procedure called node reduction ( Figure 5). Next, we show various subgraphs of the graph ofĀ to better illustrate its structure ( Figure 6). The total precedence graph ofĀ is thus the combination of Figures 5 and 6. The process of constructing the graph ofĀ starts by grouping all nodes of an event associated with a group of legs i into a single node. This can be accomplished since event nodes from the same group of legs i have "similar" incoming and outgoing arcs. As an example, consider the first set of # 1 rows ofĀ as defined in expression (53): The precedence graph for equation (73) consists of 3 × # 1 nodes, since it involves the vectors t 1 , t m , and l 1 . The relation between t 1 (k) and t 1 (k − 1) results in # 1 self connected arcs in the t 1 events with weights τ γ . Instead of expressing all elements of t 1 as individual nodes with self arcs, we reduce then to a single node with one self arc, as seen in Figure 5-a2. The dashed attribute used on the self arc indicates that for each node in the group only self arcs exist, as expressed by the "connecting" matrix E 1 . The relation between t 1 (k) and t m (k − 1) is somewhat more involved, since it contains # 1 × # m arcs, as expressed by the connecting matrix 1 1,m . The resulting node reduction is illustrated in Figure   5-b1. The node reduction for the relation between t 1 and l 1 is illustrated in Figure 5 Given the node reduction one can now proceed to construct the precedence graph ofĀ: • Figure 6-a is the graph of the block diagonal ofĀ together with the block diagonals of the sub- • Figure 6-b is the graph of the columns formed by the matrices A 12 and [A T 32 A T 42 ] T using node reductions presented in Figures 5-b1 and 5-b2.
• Figures 6-c1 and 6-c2 illustrate two subgraphs of the remaining columns ofĀ. Note that we only present the subgraphs of the first sets of # 1 and # 2 out of a total of m − 1 columns. These follow the same pattern. We use different attributes on the arcs, such as dashed, thick solid, etc., to distinguish the different node reductions, as presented in Figures 5-c1 to 5-c4.
Consider the following assumption: If assumption A2 is verified then the critical graph of G c (A) (and G c (Ā)) has a single strongly connected subgraph.
Proof: We consider two cases: In this situation the circuits presented in Figures 5-a1 and 5-a2 all belong to the critical graph since their weights are τ γ or τ ⊗m δ both equal to the max-plus eigenvalue λ. Note that any circuit c 1 of length l made from the nodes of t m , illustrated in Figure 5-a1, has an average weight of and as such also belongs to the critical graph. Any other circuit in the precedence graph ofĀ must pass through at least one node of t m , as illustrated in Figures 6-b, 6-c1, and 6-c2 (with the exception of the self-loops in Figure 5-a3 and the circuits in Figure 5-a4 that we don't consider since their weights are e and τ γ /2 both less then λ). Additionally, arcs starting in nodes from a group t q with q < m are only connected to nodes in t q+p for p 0 (or l q+p ). This is again illustrated in Figures 6-a, 6-c1, and 6-c2.
Let t [ q ] i denote element i of t q . Consider the circuit · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · a1) b1) with q < m. The average weight is (with τ γ = τ ⊗m Circuit c 2 is thus also in the critical graph. For the general circuit of the type with q 1 < q 2 < · · · < q l < m, the average weight is Again, circuit c 3 is part of the critical graph. Any circuit that passes through any node in l q , for any q, will never be in the critical graph. This is due to the fact that arcs within touchdown nodes of different leg groups yield a higher weight: As such, a path that connects a touchdown node to a lift off node "loses" τ γ − τ g = τ f from the maximum possible weight, a path from lift off to lift off nodes loses τ γ , and a path from lift off nodes to touchdown nodes loses τ g in weight. This can also be observed in the structure ofĀ, in equation (51), where the sub-matrix τ f ⊗ (τ g ⊗ W ⊕ V ) overcomes the sub-matrices τ g ⊗ W ⊕ V , τ f ⊗ W , and W . Consider, for example, the circuit c 4 : then Since all the nodes in the critical graph are connected (they are all touchdown nodes) we conclude that for the case τ γ = τ ⊗m δ = λ the critical graph ofĀ has a single strongly connected subgraph. Figure 7-a illustrates the complete critical graph ofĀ for this case.
In this situation only circuits of the type c 1 are part of the critical graph. Circuits of the type c 2 or c 3 are not part of the critical graph. Figure 7-b illustrates the resulting critical graph ofĀ. Since all the nodes of t m are connected to each other we conclude that for the case τ γ < τ ⊗m δ = λ the critical graph ofĀ has a single strongly connected subgraph.
A third case can be considered: τ ⊗m δ < τ γ = λ. In this situation the critical graph ofĀ does not have a single strongly connected subgraph. Figure 7-c illustrates this situation, that we document here without proof.
Theorem 15. Given assumptions A1 and A2, the max-plus eigenvalue λ of the system matrix A (andĀ), defined by equation (54), is unique, and the max-plus eigenvector v of A (andĀ), defined by equations (55)-(56) is unique up to a max-plus scaling factor.
Proof: According to Lemma 12 A is irreducible, and as such it has a unique max-plus eigenvalue.
According to Lemma 14 the critical graph of G c (A) has a single strongly connected subgraph, and as such its max-plus eigenvector is unique up to a max-plus scaling factor (see [1], Theorem 3.101).

E. Coupling time
Theorem 6 on page 9 describes an important property of max-plus-linear systems when the system matrix A is irreducible: it guarantees the existence of an autonomous steady-state regime that is achieved in a number of finite steps k 0 , called the coupling time. Computing the coupling time is very important for the application of legged locomotion since it provides the number of steps a robot needs to take to reach steady state after a gait transition or a perturbation.
Proof: Computing successive products ofĀ and taking advantage of its structure and equations (43) Since λ τ γ it is sufficient to show that This can be confirmed by inspecting equations (48) and (92): 1) All the terms in the upper block triangle of τ γ ⊗ W are ε while for τ f ⊗ V ⊗ W they are positive numbers.

3) In the lower block triangle
thus concluding that the coupling time is k 0 = 2 with cyclicity c = 1.

V. CONCLUSIONS
We have shown that max-plus linear models are very well suited to model leg synchronization in legged locomotion. The abstraction of the continuous-time dynamics associated with each leg as a twoevent circuit simplifies the construction of a gait space. The combinatorial nature of this gait space, arising from all the possible arrangements in which multiple legs can be synchronized, is captured by a compact representation as an ordered set of ordered sets.
For the classes of switching max-plus linear systems developed in this paper important structural properties of the max-plus system matrix are obtained in closed-form. The unique max-plus eigenvalue represents the total cycle time, the unique (up to a scaling factor) max-plus eigenvector dictates a unique steady-state behavior, and the coupling time reveals the transient response to gait switching or disturbances.
Although some of the derivations in this paper are quite lengthy, using a combination of discreteevent algebraic tools with graph-theoretic concepts has lead to an important result in robotics: since the coupling time was found to be two, we have shown that legged robots can switch gaits (or rhythms) or recover from large disturbances in at least two steps. In the derivation process we have found that similarity transformations can facilitate the algebraic manipulations by exposing the structure of the system matrices.
This was important to find closed-form expressions to the eigen-structure of the system and the coupling time, that typically need to be computed via simulations of using numerical procedures. On the graph side, the node-reduction procedure has allowed depicting graphs that can have an arbitrary large number of nodes. A graph-theoretic proof is needed for proving the uniqueness of the max-plus eigenvector. Our results are valid for robots with an arbitrary number of legs.
Further research will look towards relaxing the structure of the system matrix to address the synchronization of general cyclic systems and towards the modeling of more general gaits.