Analysis and control of max-plus linear discrete-event systems: An introduction

The objective of this paper is to provide a concise introduction to the max-plus algebra and to max-plus linear discrete-event systems. We present the basic concepts of the max-plus algebra and explain how it can be used to model a specific class of discrete-event systems with synchronization but no concurrency. Such systems are called max-plus linear discrete-event systems because they can be described by a model that is “linear” in the max-plus algebra. We discuss some key properties of the max-plus algebra and indicate how these properties can be used to analyze the behavior of max-plus linear discrete-event systems. Next, some control approaches for max-plus linear discrete-event systems, including residuation-based control and model predictive control, are presented briefly. Finally, we discuss some extensions of the max-plus algebra and of max-plus linear systems.


Introduction
In recent years both industry and the academic world have become more and more interested in techniques to model, analyze, and control complex discrete-event systems (DESs) such as flexible manufacturing systems, telecommunication networks, multiprocessor operating systems, railway networks, traffic control systems, logistic systems, intelligent transportation systems, computer networks, multi-level monitoring and control systems, and so on.Although in general DESs lead to a nonlinear description in conventional algebra, there exists a subclass of DESs for which this model becomes "linear" when it is formulated in the max-plus algebra (Baccelli et al. 1992;Cuninghame-Green 1979;Heidergott et al. 2006;Butkovič 2010), which has maximization and addition as its basic operations.More specifically, DESs 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).This leads to a description that is "linear" in the max-plus algebra.Therefore, DESs with synchronization but no concurrency are called max-plus linear DESs.
In the early sixties the fact that certain classes of DESs can be described by maxlinear models was discovered independently by a number of researchers, among whom Cuninghame-Green (1961) and Cuninghame-Green (1962) and Giffler (1960), Giffler (1963), and Giffler (1968).An account of the pioneering work of Cuninghame-Green on max-plus-algebraic system theory for DESs has been given in (Cuninghame-Green 1979).Related work on dioid theory and its applications has been undertaken by Gondran and Minoux (1976), Gondran and Minoux (1984b), and Gondran and Minoux (1987).In the late eighties and early nineties the topic attracted new interest due to the research of Cohen et al. (1985), Cohen et al. (1989), Olsder (1986), Olsder and Roos (1988), and Olsder et al. (1990a), and Gaubert (1990Gaubert ( , 1992Gaubert ( , 1993)), which resulted in the publication of Baccelli et al. (1992).Since then, several other researchers have entered the field.
The class of DESs that can be described by a max-plus linear time-invariant model is only a small subclass of the class of all DESs.However, for max-plus linear DESs there are many efficient analytic methods available to assess the characteristics and the performance of the system since one can use the properties of the max-plus algebra to analyze maxplus linear models in a very efficient way (as opposed to, e.g., computer simulation where, before being able to determine the steady-state behavior of a given DES, one may first have to simulate the transient behavior, which in some cases might require a large amount of computation time).
As will be illustrated later on in the paper, there exists a remarkable analogy between the basic operations of the max-plus algebra (maximization and addition) on the one hand, and the basic operations of conventional algebra (addition and multiplication) on the other hand.As a consequence, many concepts and properties of conventional algebra also have a max-plus analogue.This analogy also allows to translate many concepts, properties, and techniques from conventional linear system theory to system theory for max-plus linear DESs.However, there are also some major differences that prevent a straightforward translation of properties, concepts, and algorithms from conventional linear algebra and linear system theory to max-plus algebra and max-plus linear system theory for DESs.Hence, there is a need for a dedicated theory and dedicated methods for max-plus linear DESs.
In this paper we give an introduction to the max-plus algebra and to max-plus linear systems.We highlight the most important properties and analysis methods of the max-plus algebra, discuss some important characteristics of max-plus linear DES, and give a concise overview of performance analysis and control methods for max-plus linear DESs.More extensive overviews of the max-plus algebra and max-plus linear systems can be found in Baccelli et al. (1992), Cuninghame-Green (1979), Gaubert (1992), Heidergott et al. (2006), and Hardouin et al. (2018).The history of how max-plus algebra became an important tool in discrete event systems since the late seventies is described in Komenda et al. (2018).
The main feature of the current survey compared to these previous works is its compactness and its focus on analysis and model-based control for max-plus linear systems, in particular residuation-based control and model predictive control.We also include an extensive qualitative comparison between residuation-based control and model predictive control for max-plus linear systems.In addition, we provide several worked examples for basic max-plus concepts, we include several references to recent literature, and we present some results not included in previous surveys (such as, e.g., two-sided systems of linear maxplus equations, systems of max-plus-algebraic polynomial equations and inequalities, and model-based predictive control for max-plus linear systems).

Basic operations of the max-plus algebra
The basic operations of the max-plus algebra (Baccelli et al. 1992;Cuninghame-Green 1979;Heidergott et al. 2006) 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 ε def = R∪{−∞}.The reason for using these symbols is that there is a remarkable analogy between ⊕ and conventional addition, and between ⊗ and conventional multiplication: many concepts and properties from linear algebra (such as the Cayley-Hamilton theorem, eigenvectors and eigenvalues, and Cramer's rule) can be translated to the maxplus algebra by replacing + by ⊕ and × by ⊗ (see, e.g., Baccelli et al. (1992, Chapters 2, 3); Heidergott et al. (2006;Chapters 2, 5) Cuninghame-Green (1979), Gaubert (1992), and Olsder and Roos (1988)).Therefore, we also call ⊕ the max-plus-algebraic addition, and ⊗ the max-plus-algebraic multiplication.Note however that one of the major differences between conventional algebra and max-plus algebra is that in general there do not exist inverse elements with respect to ⊕ in R ε .The zero element for ⊕ is In the sequel we denote the set of non-negative integers by N = {0, 1, 2, . ..}.Let r ∈ R. The rth max-plus-algebraic power of x ∈ R is denoted by x ⊗ r and corresponds to rx in conventional algebra.If x ∈ R then x ⊗ 0 = 0 and the inverse element of x w.r.t.⊗ is x ⊗ −1 = −x.There is no inverse element for ε w.r.t.⊗ 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-algebraic operators correspond to those of conventional algebra.So max-plus-algebraic power has the highest priority, and max-plus-algebraic multiplication has a higher priority than max-plus-algebraic addition.

Max-plus-algebraic matrix operations
The basic max-plus-algebraic operations are extended to matrices as follows. (1) for all i, j .Note the analogy between (1)-( 2) and the definitions of matrix sum and product in conventional linear algebra.2), we have: The matrix ε m×n is the m × n max-plus-algebraic zero matrix: (ε m×n ) ij = ε for all i, j ; and the matrix E n is the n × n max-plus-algebraic identity matrix: (E n ) ii = 0 for all i and (E n ) ij = ε for all i, j with i = j .If the size of the max-plus-algebraic identity matrix or the max-plus-algebraic zero matrix is not specified, it should be clear from the context.The max-plus-algebraic matrix power of A ∈ R n×n ε is defined as follows: A ⊗ 0 = E n and Olsder and Roos (1988) have introduced a link between conventional algebra and the maxplus algebra based on asymptotic equivalences that allows to establish an analogy between the basic operations of the max-plus algebra (max and +) on the one hand, and the basic operations of conventional algebra (addition and multiplication) on the other hand.As a result, many concepts and properties of conventional algebra also have a max-plus analogue.In particular, Olsder and Roos (1988) used this link to show that every matrix has at least one max-plus-algebraic eigenvalue and to prove a max-plus-algebraic version of Cramer's rule and of the Cayley-Hamilton theorem.In addition, this analogy allows to translate many concepts, properties, and techniques from conventional linear system theory to system theory for max-plus linear DESs.

Connection with conventional algebra via exponentials
In De Schutter and De Moor (1997) the link introduced by Olsder and Roos (1988) has been extended and formalized.Now we recapitulate the reasoning of De Schutter and De Moor (1997) but in a slightly different form that is mathematically more rigorous.
First we extend the conventional definition of asymptotic equivalence such that we can also allow asymptotically equivalence to 0. Recall that f is asymptotically equivalent to g in the neighborhood of ∞, denoted by f (s) ∼ g(s), s → ∞, if lim s→∞ f (s) g(s) = 1.This definition in principle requires that there is no real number K such that g is identically zero in [K, ∞).However, we also say a function f is asymptotically equivalent to 0 in the neighborhood of ∞: f (s) ∼ 0, s → ∞ if there exists a real number L such that f (s) = 0 for all s L.
In this section we consider exponentials of the form e νs with s > 0. Since we want to allow exponents that are equal to ε, we set e εs equal to 0 for all positive real values of s by definition.For all x, y, z ∈ R ε we now have x ⊕ y = z ⇔ e xs + e ys ∼ (1 + δ x=y )e zs , s → ∞ (3) x ⊗ y = z ⇔ e xs • e ys = e zs for all s > 0 ( 4 ) where δ x=y = 0 if x = y and δ x=y = 1 if x = y.The relations (3) and (4) show that there exists a connection between the operations ⊕ and ⊗ performed on elements of R ε and the operations + and × performed on exponentials.
Definition 1 (Precedence graph) Consider A ∈ R n×n ε .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 a ij for each a ij = ε.
It easy to verify that every weighted directed graph corresponds to the precedence graph of an appropriately defined matrix with entries in R ε .
Now we can give a graph-theoretic interpretation of the max-plus-algebraic matrix power.Let A ∈ R n×n ε .If k ∈ N \ {0} then we have for all i, j .Hence, (A ⊗ k ) ij is the maximal weight1 of all paths of G(A) of length k that have j as their initial vertex and i as their final vertex -where we assume that if there does not exist a path of length k from j to i, then the maximal weight is equal to ε by definition.
Example 2 Consider matrix A defined in Example 1.The precedence graph G(A) of A is given in Fig. 1.Let k = 2.By direct computation (cf.Example 1), we get Now we can check that (A ⊗ 2 ) ij is the maximal weight of all paths of G(A) of length 2 that have j as their initial vertex and i as their final vertex.These paths and their corresponding weights are shown in Table 1.As one can see, the maximum weights are equal to the entries of A ⊗ 2 .
A directed graph G is called strongly connected if for any two different vertices i, j of the graph, there exists a path from i to j .

Definition 2 (Irreducible matrix)
If we reformulate this in the max-plus algebra then a matrix ) ij = ε for all i, j with i = j , since this condition means that for two arbitrary vertices i and j of G(A) with i = j there exists at least one path (of length 1, 2, . . .or n − 1) from j to i.
Example 3 Let A be defined as in Example 1.The precedence graph G(A) of A is given in Fig. 1.Clearly, G(A) is strongly connected as there exists a path from any node in G(A) to any other node, and hence A is irreducible.
3 Some basic problems in the max-plus algebra In this section we present some basic max-plus-algebraic problems and some methods to solve them.

Max-plus-algebraic eigenvalue problem
Definition 3 (Max-plus-algebraic eigenvalue we say that λ is a max-plus-algebraic eigenvalue of A and that v is a corresponding max-plus-algebraic eigenvector of A. It can be shown that matrix A ∈ R n×n ε has at least one max-plus-algebraic eigenvalue (Baccelli et al. 1992, Section 3.2.4).However, in contrast to linear algebra, the total number (multiplicities taking into account) of max-plus-algebraic eigenvalues of an n by n matrix is in general less than n.Moreover, if a matrix is irreducible, it has only one max-plusalgebraic eigenvalue (see, e.g., Cohen et al. 1985).
The max-plus-algebraic eigenvalue has the following graph-theoretic interpretation.If λ max is the maximal average weight2 over all elementary circuits of G(A), then λ max is a max-plus-algebraic eigenvalue of A. An elementary circuit is a circuit in which no vertex appears more than once, except for the initial vertex which appears exactly twice.
There exist several efficient algorithms to determine max-plus-algebraic eigenvalues such as the Karp's algorithm (Karp 1978;Cohen et al. 1985) or the power algorithm of Cochet-Terrasson et al. (1998).
To determine the max-plus-algebraic eigenvectors corresponding to a given max-plusalgebraic eigenvalue, the following procedure can be applied (Karp 1978;Cohen et al. 1985).
First we introduce the Kleene star operator3 of the matrix A: The entries of A have the following meaning: (A ) ij is the maximal weight of any path of arbitrary length in G(A) between node j and node i.We also define λ , will be a max-plus-algebraic eigenvector of A for the eigenvalue λ.This can be verified as follows: Note that in general Example 4 Consider the (irreducible) matrix A of Examples 1 and 3.The elementary circuits of G(A) are listed in Table 2.The maximum average weight is 3. Hence, λ = 3 is a max-plus-algebraic eigenvalue of A.
Example 5 For the matrix A of Example 1 we have

It can be verified that for k 5 we have
Another generalized eigenvalue problem is considered by Cochet-Terrasson et al. (1998), who define the generalized max-plus-algebraic eigenproblem for Heidergott et al. (2006, Chapter 3) use the concept of generalized eigenmode of a regular matrix A, which is defined by the pair of vectors The vector η coincides with the cycle time vector and can be seen as an extended eigenvalue, where v still remains the eigenvector.In Subiono and van der Woude (2017) a generalized power algorithm has been presented that computes the generalized eigenmode.

Systems of max-plus linear equations
In this section we consider three types of systems of max-plus linear equations, namely

A
In general, the system of max-plus linear equations A ⊗ x = b will not always have a solution, even if A is square or if it has more columns than rows.Therefore, the concept of subsolution has been introduced (see Cuninghame-Green (1979, Chapter 14), Baccelli et al. (1992, Section 3.2.3)).

Definition 4 (Subsolution) Let
Although the system A ⊗ x = b does not always have a solution, it always possible to determine the largest subsolution if we allow components that are equal to ∞ in the solution and if we assume that ε ⊗ ∞ = ∞ ⊗ ε = ε by definition. 5The largest subsolution Note that for the largest subsolution x we have A ⊗ x b.In some cases, one may want to minimize the difference between A ⊗ x and b, i.e., to find x such that max We then have max

x
Since the operation ⊕ is not invertible, an equation of the form x = A ⊗ x ⊕ b can in general not be recast into the form Ã ⊗ x = b for some matrix Ã.
If the entries of A (see ( 5)) all belong to R ε , then the least solution of x = A ⊗ x ⊕ b is given by Baccelli et al. (1992, Section 3.2.3.1):

A
A system of two-sided max-plus linear equations can be formulated as follows Walkup and Borriello (1998): , and x ∈ R n×1 ε .Note that the ith equation in Eq. ( 8) can be expanded as The maximum solution to an arbitrary system of linear max-plus equation can be obtained using the following three steps: -translate each linear max-plus equation into a small set of upper bound constraints, each of which bounds the values of a single variable from above (see Walkup and Borriello 1998, Section 2.1).-employ the max-plus closure operation to find the maximum solution to a special subset of the upper bound constraints (see Walkup and Borriello 1998, Section 2.2). -use that subset's maximum solution to guide the choice of a new constraint subset which will have a smaller maximum solution (see Walkup and Borriello 1998, Section 2.3).
The last two steps are repeated until either the process converges upon a solution which meets all the upper bound constraints, or it is found that the systems is infeasible since some variable has a maximum solution of ε.
The specific case A⊗x = C⊗x has been considered in Cuninghame-Green and Butkovic (2003).

Systems of max-plus-algebraic multivariate polynomial equations and inequalities
A system of multivariate polynomial equations and inequalities in the max-plus algebra is defined as follows: Given a set of integers {m k } k∈K and sets of coefficients {a ki } k∈K, i∈I , {b k } k∈K and set of exponents {c kij } k∈K, i∈I, j∈J where I = {1, . . ., m k }, J = {1, . . ., n} and K = {1, . . ., p eq , p eq + 1, . . ., p eq + p ineq }, find x ∈ R n such that Note that the exponents c kij can be negative or real.In conventional algebra the above equations can be written as In De Schutter and De Moor (1996), De Schutter (1996), and De Schutter and De Moor (1998) it has been shown that the above problem and related max-plus problems such as computing max-plus matrix decompositions, transformation of max-plus linear state space models, state space realization of max-plus linear systems, construction of matrices with a given max-plus characteristic polynomial, and solving systems of max-min-plus equations can be recast as a so-called extended linear complementarity problem (ELCP), which is defined as follows:  DESs with only synchronization and no concurrency can be modeled by a max-plusalgebraic model of the following form Baccelli et al. (1992), Cuninghame-Green (1979), and Heidergott et al. (2006): , and C ∈ R l×n ε , where m is the number of inputs and l the number of outputs.The vector x represents the state, u is the input vector, and y is the output vector of the system.It is important to note that in ( 10)-( 11) the components of the input, the output, and the state are event times, and that the counter k in ( 10)-( 11) is an event counter.For a manufacturing system (see also Example 7 below), u(k) typically represents the time instants at which raw material is fed to the system for the kth time, x(k) the time instants at which the machines start processing the kth batch of intermediate products, and y(k) the time instants at which the kth batch of finished products leaves the system.
Due to the analogy with conventional linear time-invariant systems, a DES that can be modeled by ( 10)-( 11) will be called a max-plus linear time-invariant DES.
Typical examples of systems that can be modeled as max-plus linear DESs are production systems, railroad networks, urban traffic networks, queuing systems, and legged robots (Baccelli et al. 1992;Cuninghame-Green 1979;Heidergott et al. 2006;Lopes et al. 2014).We will now illustrate in detail how the behavior of a simple manufacturing system can be described by a max-plus linear model of the form ( 10)-( 11).
Example 7 Consider the system of Fig. 2.This production system consists of three processing units: P 1 , P 2 , and P 3 .Raw material is fed to P 1 and P 2 , processed, and sent to P 3 where assembly takes place.The processing times for P 1 , P 2 , and P 3 are respectively d 1 = 12, d 2 = 11, and d 3 = 7 time units.We assume that it takes t 2 = 2 time units for the raw material to get from the input source to P 2 and that it takes t 4 = 1 time unit for the finished products of processing unit P 2 to reach P 3 .The other transportation times (t 1 , t 3 , and t 5 ) are assumed to be negligible.We assume that at the input of the system and between the processing units there are buffers with a capacity that is large enough to ensure that no buffer overflow will occur.We consider the situation where initially all buffers are empty and none of the processing units contains raw material or intermediate products.This corresponds to in fact to the case where x(0) = ε 3×1 ; if initially, some buffers or some processing units are non-empty, then we will have x(0) = ε 3×1 .
A processing unit can only start working on a new product if it has finished processing the previous product.We assume that each processing unit starts working as soon as all parts are available.Define u(k): time instant at which raw material is fed to the system for the kth time, x i (k): time instant at which the ith processing unit starts working for the kth time, y(k): time instant at which the kth finished product leaves the system.
Let us now determine the time instant at which processing unit P 1 starts working for the kth time.If we feed raw material to the system for the kth time, then this raw material is available at the input of processing unit P 1 at time t = u(k) + 0. However, P 1 can only start working on the new batch of raw material as soon as it has finished processing the previous, i.e., the (k − 1)st, batch.Since the processing time on P 1 is d 1 = 12 time units, the (k − 1)st Discrete Event Dynamic Systems (2020) 30:25-54 intermediate product will leave P 1 at time t = x 1 (k − 1) + 12. Since P 1 starts working on a batch of raw material as soon as the raw material is available and the current batch has left the processing unit, this implies that we have for k = 1, 2, . . .The condition that initially processing unit P 1 is empty and idle corresponds to the initial condition x 1 (0) = ε since then it follows from (12) that x 1 (1) = u(1), i.e., the first batch of raw material that is fed to the system will be processed immediately.
Using a similar reasoning we find the following expressions for the time instants at which P 2 and P 3 start working for the (k + 1)st time and for the time instant at which the kth finished product leaves the system: for k = 1, 2, . . .The condition that initially all buffers are empty corresponds to the initial condition Let us now rewrite the evolution equations of the production system using the symbols ⊕ and ⊗.It is easy to verify that (12) can be rewritten as If we also do this for ( 13)-( 15) and if we rewrite the resulting equations in max-plusalgebraic matrix notation, we obtain T .Note that this is a model of the form (10)-( 11).
In the next section we shall use this production system to illustrate some of the maxplus-algebraic techniques that can be used to analyze max-plus linear time-invariant DESs.

Analysis of max-plus linear systems
Now we present some analysis techniques for max-plus linear DESs that can be described by a model of the form (10)-( 11).
First we determine the input-output behavior of the max-plus linear DES.We have Discrete Event Dynamic Systems (2020) 30:25-54 etc., which yields k=1 be the output sequence that corresponds to the input sequence u 1 (with initial condition x 1 (0) = x 1,0 ) and let y 2 = {y 2 (k)} ∞ k=1 be the output sequence that corresponds to the input sequence u 2 (with initial condition x 2 (0) = x 2,0 ).Let α, β ∈ R ε .From ( 16) it follows that the output sequence that corresponds to the input sequence ) is given by α ⊗ y 1 ⊕ β ⊗ y 2 .This explains why DESs that can be described by a model of the form (10)-( 11) are called max-plus linear.
Let p ∈ N \ {0}.If we define with For the production system of Example 7 this means that if we know the time instants at which raw material is fed to the system and we know the initial state x(0), we can compute the time instants at which the finished products will leave the system.Often we assume that x(0) = ε n×1 .For the simple production system of Example 7 this would mean that all the buffers are initially empty.
Example 8 Consider the production system of Example 7. Define Y = y(1) y( 2

) y(3) y(4)
T and U = u(1) u( 2) u(3) u( 4) If we feed raw material to the system at time instants u(1) = 1, u(2) = 8, u(3) = 15, u(4) = 19, the finished products will leave the system at time instants y(1) = 22, y(2) = 33, y(3) = 44, and y(4) = 56 since Now we consider the autonomous max-plus linear DES described by with x(0) = x 0 .For the production system of Example 7 this would mean that we start from a situation in which some internal buffers and all the input buffer are not initially empty (if x 0 = ε n×1 ) and that afterwards the raw material is fed to the system at such a rate that the input buffers never become empty.
If the system matrix A of the autonomous DES is irreducible, then we can compute the "ultimate" behavior of the autonomous DES by solving the max-plus-algebraic eigenvalue problem for i = 1, 2, . . ., n and for all k k 0 .From this relation it follows that for a production system the average time duration of a cycle of the production process when the system has reached its cyclic behavior will be given by λ.The average production rate will then be equal to 1 λ .This also enables us to calculate the utilization levels of the various machines in the production process.Furthermore, some algorithms to compute the eigenvalue also yield the critical paths of the production process and the bottleneck machines (Cohen et al. 1985).

Example 9
The system matrix A of the production system of Example 7 is not irreducible and it does not lead to a behavior of the form 6 (18).In fact, it can be verified that A has three eigenvalues λ 1 = 12, λ 2 = 11, and λ 3 = 7, with corresponding eigenvectors So by choosing k 0 = 1, we see that for k ≥ k 0 the first and the third element of x satisfies x i (k + c) − x i (k) = cλ 1 for i = 1 or 3 with c = 1, but the second element of x satisfies x 2 (k + c) − x 2 (k) = cλ 2 with c = 1, which means that (18) does not hold in this case. 6Due to the structure of the system matrix A it is easy to verify that A ⊗ k 11 12 ⊗ k and A ⊗ k 22 = 11 ⊗ k for all k 1.Hence, the relation given in Theorem 1 does not hold in this case.

Control of max-plus linear DES
The basic control problem for max-plus linear DESs consists in determining the optimal input times (e.g., feeding times of raw material or starting times of processes or activities) for a given reference signal (e.g., due dates for the finished products or completion dates for processes or activities).In the literature many different approaches are described to solve this problem.Among these the most common ones are based on residuation and on model predictive control (MPC).Residuation essentially consists in finding the largest solution to a system of max-plus inequalities with the input times as variables and the due dates as upper bounds.The MPC approach is essentially based on the minimization of the error between the actual output times and the due dates, possibly subject to additional constraints on the inputs and the outputs.
Remark 1 For simplicity, we only consider single-input single-output (SISO) systems in this section.Note however MPC can very easily be extended to multi-input multi-output (MIMO) systems.

Residuation-based control
The basic control problem for max-plus linear DESs consists in determining the optimal feeding times of raw material to the system and/or the optimal starting times of the (internal) processes.
Consider ( 17) with x(0) = ε n×1 .If we know the vector Y of latest times at which the finished products have to leave the system, we can compute the vector U of (latest) time instants at which raw material has to be fed to the system by solving the system of max-plus linear equations H ⊗ U = Y , if this system has a solution, or by determining the largest subsolution of H ⊗ U = Y , i.e., determining the largest U such that H ⊗ U Y .This approach is also based on residuation (Blyth and Janowitz 1972).
Note that the method above corresponds to just-in-time production.However, if we have perishable goods it is sometimes better to minimize the maximal deviation between the desired and the actual finishing times.In that case we have to solve the problem min This problem can be solved using formula (7).Note that the largest deviation between the desired finishing and the actual finishing times is now equal to δ 2 = 2, which means that the maximal deviation between the desired (Y ) and the actual ( Ỹ ) finishing times is minimized.
In particular, Libeaut and Loiseau (1995) have applied the geometric approach and residuation theory in order to find the optimal input.The geometric approach is a collection of mathematical concepts developed to achieve a better and neater insight into the most important features of multi-variable linear dynamical systems in connection with compensator and regulator synthesis problems.It is based on the state space representation and it also easily links SISO and MIMO systems and clarifies in a concise and elegant way some common properties that cannot be obtained by the transform-based techniques usually adopted in the SISO case.Related work can be found in Ilchmann (1989).Using these results, in Libeaut and Loiseau (1995) the set of admissible initial conditions of a linear system is defined and characterized geometrically and the optimal input is computed by applying residuation theory.In Boimond and Ferrier (1996) the Internal Model Control (IMC) structure used in conventional control theory is extended to deterministic max-plus linear DESs.The IMC philosophy relies on the internal model principle, which states that control can be achieved only if the control system encapsulates, either implicitly or explicitly, some representation of the process to be controlled; a comprehensive explanation can be found in Garcia and Morari (1982).The controller design raises the problem of model inversion, where the residuation approach also plays an important role.In Menguy et al. (1997), a feedback control structure is proposed to be able to take into account a possible mismatch between the system and its model.Instead of adopting the commonly used IMC approach for closed-loop systems, the authors proposed another closed-loop control structure consisting in applying an open-loop control approach that is modified by using the system output behavior.In fact, the model is initially supposed to be exact; subsequently, the control structure is modified by using the system output in order to be as close as possible to the optimal system control.The optimization problem is solved using residuation.The method of Menguy et al. (1998) is also based on inverting a dynamic system and applying residuation theory.The proposed control structure is based on adaptive control and encompasses both identification and inversion of a dynamic system.In Lahaye et al. (2008), a just-in-time optimal control for a first-in first-out (FIFO) time event graph is proposed based on residuation theory.The aim is to compute the largest control u such that the firing dates of output transitions (or simply the output signals) occur at the latest before the desired ones.In Brunsch et al. (2012), Brunsch andRaisch (2012), andDavid-Henriet et al. (2017) residuation-based control is applied for max-plus linear systems arising in the context of manufacturing systems and of high-throughput screening (e.g., for the pharmaceutical industry).
Note that the sequence u(1), u(2), . . ., u(N) of the feeding times should be nondecreasing as it corresponds to a sequence of consecutive feeding times.However, in general a residuation-based solution does not always satisfy this property.This problem can be overcome by projection onto the set of non-increasing sequences as in the approach proposed by Menguy et al. (2000a).In addition, an initial state x(0) that is not necessarily equal to ε n×1 can be included and an explicit form can be derived for the residuation controller.The resulting expression involves the operations minimization and addition from the minplus algebra, which is the dual of the max-plus algebra and which has minimization (⊕ ) and addition (⊗ ) as basic operations with the zero element ε = ∞.More specifically, the residuation controller is now given by the following expression (van den Boom and De Schutter 2014): where U 0 and S are respectively a vector with p components and a p by p matrix defined as7

Model predictive control
A somewhat more advanced control approach for max-plus linear DESs has been developed by De Schutter and van den Boom (2001).This approach is an extension to max-plus linear DESs of the model-based predictive control approach called Model Predictive Control (MPC) (Camacho and Bordons 1995;Maciejowski 2002;Rawlings and Mayne 2009) that has originally been developed for time-driven systems.
The main advantage of the MPC method of De Schutter and van den Boom ( 2001), comparing to other available methods for control design in max-plus linear DES, is that it allows to include general linear inequality constraints on the inputs, states, and outputs of the system.
In MPC for max-plus linear DESs at each event step k the controller computes the input sequence u(k), . . ., u(k + N p − 1) that optimizes a performance criterion J over the next N p event steps, where N p is called the prediction horizon, subject to various constraints on the inputs, states, and outputs of the system.Typically, the performance criterion J = J out + λJ in aims at minimizing the difference or the tardiness with respect to a due date signal r(•), while at the same time making the inputs as large as possible (just-in-time production) where λ ∈ (0, 1) is a weighting factor.A typical choice for the output cost function is and for the input criterion function a typical choice is Note that J out penalizes late delivery and J in penalizes early feeding times (i.e., it favors just-in-time feeding).
This results in an optimization problem that has to be solved at each event step k.In order to reduce the computational complexity, often a control horizon N c is introduced with N c < N p and it is assumed that the feeding rate ( u(k + j ) = u(k + j ) − u(k + j − 1)) is constant after event step k+N c , i.e., u(k+N c + ) = u(k+N c −1) for = 0, . . ., N p −N c −1.Moreover, we can also have lower and upper bound constraints on x(k + j ), y(k + j ), y(k + j ), or linear constraints on x(k + j ), y(k + j ), and u(k + j ).In general, this leads to a system of linear inequalities Hence, at event step k the following optimization problem has to be solved: 24) MPC uses a receding horizon approach.This means that once the optimal input sequence has been determined only the input for the first event step is applied to the system.Next, at event step k + 1 the new state of the system is determined or measured8 , the prediction window is shifted, and the whole process is repeated again.This receding horizon approach introduces a feedback mechanism, which allows to reduce the effects of possible disturbances and model mismatch errors.More information on MPC for max-plus-linear systems can be found in De Schutter and van den Boom (2001); Goto (2009); Necoara et al. (2007Necoara et al. ( , 2009a)); van den Boom and De Schutter (2002a).MPC for max-plus-linear systems with partial synchronization is proposed in (David-Henriet et al. 2016).

Comparison of residuation-based control and MPC for max-plus-linear systems
In this section we compare the control methods based on residuation with the ones based on MPC.We discuss four items: constraint handling, cost functions, computation time, and implementation.We end the section with a worked example.

Constraint handling
First note that the problem of solving the problem ( 23)-( 25) with output criterion (21) and input criterion ( 22), but without constraints ( 26)-( 28), is equivalent to solving the residuation problem with p = N p in a receding horizon setting.If we also take constraint (26) into account, we obtain the result from Menguy et al. (2000a) with a non-decreasing input signal, for which the explicit solution is given by (20).
As was already mentioned in the previous section, the main advantage of the MPC, compared to other available methods for control design for max-plus linear DES, is that it allows to include general linear inequality constraints on the inputs, states, and outputs of the system.This is done by writing the control law as the result of a constrained optimization problem ( 23)-( 28).MPC is until so far the only known method to handle the constraint (28).

Different cost functions
In Section 5.2.2 we have discussed MPC with the performance criterion J = J out + λJ in where J out and J in are given by ( 21) and ( 22), respectively.This performance criterion results in a just-in-time control strategy, which is comparable to the ones in residuationbased control.
However, in (De Schutter and van den Boom 2001) it has been shown that MPC can handle a broad range of other performance criteria.Some alternative output criteria are where (29) does not only penalize late delivery but also any deviation from the due date, and (30) can be used if we want to balance the output rates.Some alternative input criteria are where (31) minimizes the maximum holding time of products in the system, and (32) balances the input rates.

Computation time
An important advantage of residuation-based control w.r.t.MPC is that of the constraint (28) is not present residuation-based control offers an analytic solution that can be computed very efficiently, while in MPC at every event step an optimization problem has to be solved numerically.
In some cases (De Schutter and van den Boom 2001) the MPC optimization will result in a linear programming problem, which can to be solved numerically and on-line in a finite number of iterative steps using reliable and fast algorithms.This is the case if we, e.g., solve the problem ( 23)-( 28) with output criterion (21) and input criterion ( 22), and if the matrix B c in constraint (28) satisfies [B c ] ij ≥ 0 for all i, j .
If the performance criterion is not convex the problem has to be solved with other optimization techniques, such as an ELCP algorithm (De Schutter and De Moor 1995), a mixed integer linear programming problem (Heemels et al. 2001;Bemporad and Morari 1999), or a optimistic programming algorithm (Xu et al. 2014).
The time required for the optimization makes model predictive control not always suitable for fast systems and/or complex problems.

Implementation
Another difference between residuation-based controllers and controllers based on MPC is the complexity in implementation.
The control law of a residuation-based controller can be written as an analytic expression, which can be implemented on a programmable logic controller (PLC), a distributed control system (DCS), or a programmable automation controller (PAC) in a straightforward way.
Most MPC controllers in industry use personal computers or dedicated microprocessorbased systems to manipulate the data and to perform the calculations, usually needing dedicated optimization software.This means that the final step in the controller design (implementation) is more complex for MPC controllers than for residuation-based controllers.
For application in industry it is important that controllers are cheap and therefore the additional performance and scope of the MPC controller should weigh up against the higher implementation costs.
The residuation-2 method (Libeaut and Loiseau 1995;Menguy et al. 1997) includes the non-decreasing input constraint and it yields a non-decreasing -and thus feasible -input sequence.Note that the residuation-2 approach is equivalent to solving MPC problem ( 23)-(26) for N p = 15.
For k ∈ {8, 9, 10} we see observe that the tracking error in the MPC approach is negative.This due to the fact that the increment input signal hits the constraint u(k) ≤ 15.The residuation methods both show input signals with an increment larger than 15, which thus do not satisfy all the imposed constraints.
6 Related work in modeling, performance analysis, identification, and control Addad et al. (2010) present an approach to evaluate the response time in networked automation systems that use a client/server protocol.The developments introduced are derived from modeling the entire architecture in the form of timed event graphs, as well as from the resulting state space representation in max-plus algebra.Another method for deriving a max-plus linear state space representation for repetitive FIFO systems is presented in Goto and Masuda (2008).

Modeling and performance analysis
An interesting topic is the use of network calculus as a tool to analyze the performance in communication networks and queuing networks, in particular to obtain deterministic bounds.Although network calculus is originally based on min-plus algebra, alternative formulations can be developed based on max-plus algebra (Liebeherr 2017).In Bouillard and Thierry (2008) some efficient algorithms, implementing network calculus operations for some classical functions, have been provided as well as the analysis of their computational complexity.
In (van den Boom and De Schutter 2006;2012) switching max-plus linear systems were introduced as an extension of max-plus linear systems.The system can switch between different modes, where in each mode the system is described by a max-plus-linear state equation and a max-plus-linear output equation, with different system matrices for each mode.

Identification and verification
In Schullerus et al. (2006) the problem of designing adequate input signals for state space identification of max-plus linear systems is addressed.This work constitutes an improvement compared to the existing methods by adding additional constraints due to safety or performance requirements on input and output signals besides reducing the computational burden of the already existing models.
Observer design for max-plus linear systems is considered in Hardouin et al. (2010) and Hardouin et al. (2017).Stochastic filtering of max-plus linear systems with bounded disturbances is discussed in Santos-Mendes et al. (2019).Adzkiya et al. (2013) develop a method to generate finite abstractions (i.e., simplified representations that still capture a given behavior or feature of the original system) of max-plus-linear models.This approach enables the study -in a computationally efficient way -of general properties of the original max-plus-linear model by verifying (via model checking) equivalent logical specifications over the finite abstraction.Related work is reported in Esmaeil Zadeh Soudjani et al. (2016).

Control
In Katz (2007) the extension of the geometric approach to linear systems over the max-plus algebra is presented.This approach is based on the state space representation rather than using residuation methods; however, it is still different from the MPC approach.The aim is to compute the set of initial states for which there exists a sequence of control vectors that makes the state of system (10) converge in the given space.Related work is described in Shang (2013) and Shang et al. (2016).
An important topic in the study of control for max-plus linear system is the incorporation of uncertainty in the system parameters.Note that noise and parameter uncertainty in max-plus linear systems will appear in a max-plus-multiplicative way as perturbations of the system parameters (Olsder et al. 1990b), usually leading to delays in the system, even if the uncertainty has a zero mean.This perturbation can have a bounded character or it can be modeled in a stochastic way.The bounded case has been studied in Lahaye et al. (1999), Lhommeau et al. (2004), van den Boom and De Schutter (2002b), and Necoara et al. (2009b).In (van den Boom and De Schutter 2004) it is shown that the stochastic MPC problem can be recast into a convex optimization problem.To reduce the complexity of the stochastic MPC optimization problem Farahani et al. (2016) use the moments of a random variable to obtain approximate solution using less computation time.Xu et al. (2019) introduced chance constraints in the MPC problem for stochastic max-plus linear systems.

Summary
This paper has presented a survey of the basic notions of the max-plus algebra and maxplus linear discrete-event systems (DESs).We have introduced the basic operations of the max-plus algebra and presented some of the main definitions, theorems, and properties of the max-plus algebra.Next, we have given an introduction to max-plus linear DES, and discussed some elementary analysis and control methods for max-plus linear DESs including worked examples.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 1
Fig. 1 Precedence graph of the matrix A of Examples 1 and 2. The vertices are indicated in an encircled bold italic font, and the weights are indicated next to the arcs in a regular font 2 1 5 and k = 5, 6, . . .For given matrices A, B ∈ R n×n ε the generalized or two-sided max-plus-algebraic eigenproblem (Cuninghame-Green and Butkovič 2008; Gaubert and Sergeev 2013; Butkovič and Jones 2016) consists in finding λ ∈ R ε and a vector or equivalently, x = −(A T ⊗ (−b)) Example 6 Consider the matrix A of Example 1 and let b = 1 2 3 T .The system of equations A ⊗ x = b does not have a solution.However, the largest subsolution is given by x = −1 −2 0 T , and we have A ⊗ x = 1 0 3 T b.
9 ) subject to Ax c and Bx = d.Algorithms for solving ELCPs can be found in De Schutter and De Moor (1995) (for computing the entire solution set) and in De Schutter et al. (2002) (for finding only one solution).

Fig. 2
Fig. 2 A simple production system

Example 10
Let us again consider the production system of Example 7 and the matrix H and the vectors U and Y as defined in Example 8.If the finished parts should leave the system before time instants 21, 32, 48, and 55 and if we want to feed the raw material to the system as late as possible, then we should feed raw material to the system at time instants 0, 11, 23, 34 since the largest subsolution of H ⊗ U = Y = 21 32 48 55 T is Û = 0 11 23 34 T .The actual output times Ŷ are given by Ŷ = H ⊗ Û = 21 32 44 55 T .Note that Ŷ ≤ Y .The largest deviation δ between the desired and the actual output times is equal to 4. The input times that minimize this deviation are given by Ũ = Û ⊗ δ 2 = Û ⊗2 = 2 13 25 36 T .The corresponding output times are given by Ỹ = 23 34 46 57 T .
for properly defined matrices and vectors A c (k), B c (k), C c (k), and d c (k), and where

Fig. 3
Fig. 3 Tracking error y − r obtained from residuation and MPC Samira S. Farahani earned her B.Sc. degree in applied mathematics from the Sharif University of Technology, Iran, in 2005.She earned her M.Sc.degree in applied mathematics and Ph.D. degree in systems and control from the Delft University of Technology, The Netherlands, in 2008 and 2012, respectively.She is a postdoctoral researcher at the Engineering Systems and Services Department at the Delft University of Technology.

Table 2
The elementary circuits of the precedence graph of the matrix A of Examples 1, 3, and 4