Tutorial for the max-plus toolbox

In this section, we will look at the max_plus toolbox. It offers a way to carry out operations with max-plus matrices, compute trajectories for a timed Petri net, graphically show such trajectories in Gantt charts and analyze and draw precedence graphs. The first thing to do is to import the max_plus module from folder petritub:

from petritub.max_plus import *

Max-plus matrix operations

Consider the following max-plus matrices, with \(\varepsilon = - \infty\) and \(e = 0\):

\[\begin{split}\begin{align} A_0 = \begin{bmatrix} \varepsilon & \varepsilon & \varepsilon\\ 1 & \varepsilon & 3\\ 4 & \varepsilon & \varepsilon\\ \end{bmatrix},\quad A_1 = \begin{bmatrix} 7 & 1 & 2\\ 3 & \varepsilon & e\\ 1 & \varepsilon & 4\\ \end{bmatrix} \end{align}\end{split}\]

In our toolbox, a max-plus matrix (or scalar, i.e., a \(1\times 1\) matrix) is an instance of the Python class Matrix. It can be initialized by specifying its elements in a two-dimensional list:

A = Matrix([[eps,eps,eps],
            [  1,eps,  3],
            [  4,eps,eps]])

B = Matrix([[  e,  1,eps],
            [  3,eps,eps],
            [  1,eps,  e]])

To print the matrices, use function print(). This function displays \(\varepsilon\) as “ε” and \(e\) as “e”.

print("B = ",B)

B =
     e     1     ε
     3     ε     ε
     1     ε     e

Some computers, however, do not render “ε” correctly; to visualize \(\varepsilon\) as “-inf” and “e” as “0” you can use the __format__() method with parameter 'v' (or 'verbose') as follows:

print("B = ",format(B,'v'))
# or
#print("B = {:v}".format(B))

B =
      0      1   -inf
      3   -inf   -inf
      1   -inf      0

We can determine specific properties of these matrices using the methods of the library. As an example, we can obtain the size of a Matrix object and determine if it is reducible in the following way:

print(A.size())
print(A.is_reducible())

(3, 3)
True

Operations \(\oplus\) (element-wise max-plus addition) and \(\otimes\) (max-plus matrix multiplication) are called using + and *, respectively:

print("A + B = ",A + B)
print("A * B = ",A * B)

A + B =
     e     1     ε
     3     ε     3
     4     ε     e

A * B =
     ε     ε     ε
     4     2     3
     4     5     ε

We can also use the Kleene star operator to compute \(A^*\) by calling the star() method of class Matrix.

print("A* = ",A.star())

A* =
     e     ε     ε
     7     e     3
     4     ε     e

The Kleene star operator \(X^*\) can be computed (i.e., converges to a matrix of real numbers and \(\varepsilon\)’s) on a matrix \(X\) only if the precedence graph corresponding to \(X\), denoted \(\mathcal{G}(X)\), does not contain circuits with positive weight. When you use method star(), this is automatically detected and if the condition is not satisfied, a warning is returned.

For instance, graph \(\mathcal{G}(B)\) contains such a circuit; we can check this manually by drawing \(\mathcal{G}(B)\):

G_B = PrecedenceGraph(B)
G_B.draw()

The circuit \(1\rightarrow 2\rightarrow 1\) has a positive weight (of 4). This is confirmed by the computation of element \((1,1)\) of matrix \(B^2 = B\otimes B\), as \((B^2)_{11} = 4\). To raise a matrix to a power, in the max-plus sense, you can use **:

B_squared = B ** 2
print("B²(1,1) =", B_squared[0,0])     # print element (1,1) of B_squared. Note that Python uses a zero-based numbering

B²(1,1) = 4

If we try to compute \(B^*\), we will get a warning:

B.star()

Warning! The star operator of the matrix

     e     1     ε
     3     ε     ε
     1     ε     e
cannot be computed (there exist circuits with positive weight in its precedence graph)!

<petritub.max_plus.Matrix at 0x71529bbc6aa0>

The max_plus toolbox offers other functionalities for the analysis of precedence graphs; these are discussed in Section Precedence graphs.

Furthermore, it is possible to compute the largest eigenvalue of a max-plus matrix (i.e., the maximum mean weight of all circuits in the precedence graph associated to the matrix), as well as the set of its linearly independent eigenvectors corresponding to the maximum eigenvalue. For matrix \(B\), this is done as follows:

print("Largest eigenvalue of B: ", B.eigenvalue())
print("Linearly independent eigenvectors of B corresponding to the largest eigenvalue:", B.eigenvectors())

Largest eigenvalue of B:  2.0
Linearly independent eigenvectors of B corresponding to the largest eigenvalue:
       -1.0
          e
       -2.0

This means that the largest eigenvalue of \(B\) is 2 and that \(B\) has only one linearly independent eigenvector:

\[\begin{split}\begin{bmatrix} -1\\ e\\ -2\\ \end{bmatrix}\end{split}\]

Computing a trajectory

Autonomous system

We want to study the trajectory \(\{x(k)\}_{k=0,1,\ldots,5}\) of a max-plus linear dynamical system of the form

\[x(k+1) = A\otimes x(k),\]

where matrix \(A\) is

\[\begin{split}A = \begin{bmatrix} 1 & \varepsilon & e\\ 5 & e & 4\\ \varepsilon & 2 & \varepsilon \end{bmatrix},\end{split}\]

and \(x(0) = [e\ 5\ 8]'\).

To generate \(\{x(k)\}_{k=0,1,\ldots,5}\), we can call function trajectory as follows.

# matrix A
A = Matrix([[  1,eps,  e],
            [  5,  e,  4],
            [eps,  2,eps]])

# initial state x(0)
x_0 = [e, 5, 8]

# final k value
k_end = 5

# compute trajectory
tr = trajectory(A, x_0, k_end)

State trajectory:
e     8     9    14    15    20
5 -> 12 -> 13 -> 18 -> 19 -> 24
8     7    14    15    20    21

Non-autonomous system

The following figure, taken from [1], shows a manufacturing system with three machines. Each of them has a processing time and a given maximum capacity. The first two machines need to pre-process raw parts, which are then needed in the third machine for the finished product.

This manufacturing system can be modeled by the timed Petri net (more specifically, the timed event graph) shown in the following figure [1]. The timing information is encoded in the holding times for some of the places.

To compute state and output trajectories, we first obtain matrices \(A_0,\ A_1,\ A_3,\ B_0,\ C_0\) of the implicit state equation for our system

\[x(k+1) = A_0 \otimes x(k+1) \oplus A_1 \otimes x(k) \oplus A_3\otimes x(k-2) \oplus B_0 \otimes u(k+1),\]

\[y(k) = C_0 \otimes x(k).\]

Here, the internal transitions (i.e., transitions which are neither input nor output transitions) are numbered with \(x_i\) with \(i = 1,...,6\), input transitions as \(u_1\), \(u_2\), and the output transition is \(y\). In the max-plus algebra, daters \(x_i(k),\ u_i(k),\ y(k)\) are defined as the (earliest) time for the \(k\)-th firing of transitions \(x_i\), \(u_i\), and \(y\), respectively. For our above system we can obtain the following max-plus matrices:

\[\begin{split}A_0 = \begin{bmatrix} \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ 2 & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & 5 & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & 1 & \varepsilon & 3 & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & 12 & \varepsilon \end{bmatrix},\quad A_1 = \begin{bmatrix} \varepsilon & e & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & e & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon \end{bmatrix},\end{split}\]

\[\begin{split}A_3 = \begin{bmatrix} \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & e\\ \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon \end{bmatrix}, \quad B_0 = \begin{bmatrix} 1 & \varepsilon \\ \varepsilon & \varepsilon \\ \varepsilon & 2 \\ \varepsilon & \varepsilon \\ \varepsilon & \varepsilon \\ \varepsilon & \varepsilon \end{bmatrix},\quad C_0 = \begin{bmatrix} \varepsilon & \varepsilon & \varepsilon & \varepsilon & \varepsilon & e \end{bmatrix}.\end{split}\]

The Matrices are defined below.

A0 = Matrix([[eps,eps,eps,eps,eps,eps],
             [  2,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,  5,eps,eps,eps],
             [eps,  1,eps,  3,eps,eps],
             [eps,eps,eps,eps, 12,eps]])
A1 = Matrix([[eps,  e,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,  e,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps]])
A3 = Matrix([[eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,eps],
             [eps,eps,eps,eps,eps,  e],
             [eps,eps,eps,eps,eps,eps]])
B0 = Matrix([[  1,eps],
             [eps,eps],
             [eps,  2],
             [eps,eps],
             [eps,eps],
             [eps,eps]])
C0 = Matrix([[eps,eps,eps,eps,eps,  e]])

Implicit dynamical equations cannot be used directly to compute a trajectory. Therefore, we need to write them in explicit form as

\[x(k+1) = \underbrace{A_0^*A_1}_{\overline{A}_1} x(k) \oplus \underbrace{A_0^*A_3}_{\overline{A}_3} x(k-2) \oplus \underbrace{A_0^*B_0}_{\overline{B}_0} u(k+1).\]

Thus, in the following we compute \(\overline{A}_1\), \(\overline{A}_3\), and \(\overline{B}_0\).

# The following line throws a warning in the case the precedence graph of A0 has a positive circuit; in that case, the explicit form cannot be obtained
A0star = A0.star()
A1bar = A0star * A1
A3bar = A0star * A3
B0bar = A0star * B0

In the current version of the toolbox (maybe more flexibility will be added in the future), a last step is needed before computing the trajectory: rewrite the dynamics as

\[\tilde{x}(k+1) = A \tilde{x}(k) \oplus B \tilde{u}(k),\]

\[y(k) = C \tilde{x}(k),\]

where

\[\begin{split}\tilde{x}(k) = \begin{bmatrix} x(k)\\ x(k-1)\\ x(k-2) \end{bmatrix},\quad A = \begin{bmatrix} \overline{A}_1 & N & \overline{A}_3\\ E & N & N\\ N & E & N \end{bmatrix},\end{split}\]

\[\begin{split}B = \begin{bmatrix} \overline{B}_0 \\ N \\ N \end{bmatrix},\quad C = \begin{bmatrix} C_0 & N & N \end{bmatrix},\end{split}\]

\(N\) is the matrix with only \(\varepsilon\)’s and \(E\) is the identity matrix in the max-plus algebra. The two matrices can be called in the toolbox using zeros(n) (or zeros(n,n)) and eye(n), respectively, where n indicates the dimension of the matrices.

Function concatenate allows to write matrices \(A\), \(B\), \(C\) in block-diagonal form. Note that the optional parameter axis lets you decide in which “direction” to perform the concatenation: * axis=0 (default) \(\rightarrow\) vertical concatenation, * axis=1 \(\rightarrow\) horizontal concatenation.

# matrix A
A = concatenate((
        concatenate((   A1bar,zeros(6),   A3bar),axis=1),
        concatenate((  eye(6),zeros(6),zeros(6)),axis=1),
        concatenate((zeros(6),  eye(6),zeros(6)),axis=1)))
# matrix B
B = concatenate((B0bar,zeros(6,2),zeros(6,2)))
# matrix C
C = concatenate((C0,zeros(1,6),zeros(1,6)),axis=1)

We can initialize the state vector \(\tilde{x}(0) = \begin{bmatrix}x(0)'\ x(-1)'\ x(-2)'\end{bmatrix}'\) arbitrarily; as usual, we impose \(\tilde{x}_i(0) = \varepsilon\) for all \(i\), which corresponds to assuming that the \(0\)-th, \(-1\)-st, and \(-2\)-nd firings of the internal transitions have never occurred (or, equivalently, that they occurred at time \(\varepsilon=-\infty\)). This is motivated by the real example: before the entrance of the first components in the machines of the manufacturing system, the machines have never been started.

As for \(u(k)\), let us suppose, for instance, that \(u_i(k) = e\) for all \(i = 1,2\), and for all \(k\in\{1,2,\dots,5\}\); this means that \(5\) firings of the internal transitions occur at time \(0\), i.e., machines M1 and M2 receive \(5\) raw parts each at time \(0\). We store signal \(\{u(k)\}_{k=1,2,\ldots,5}\) in a matrix \(U\):

\[\begin{split}U = \begin{bmatrix} u(1)& u(2) & u(3) & u(4) & u(5) \end{bmatrix}= \begin{bmatrix} u_1(1) & u_1(2) & u_1(3) & u_1(4) & u_1(5)\\ u_2(1) & u_2(2) & u_2(3) & u_2(4) & u_2(5)\\ \end{bmatrix} = \begin{bmatrix} e & e & e & e & e\\ e & e & e & e & e\\ \end{bmatrix}.\end{split}\]

We can create a matrix of dimensions \(m\times n\) populated only by \(e\)’s using function ones(m,n). The name of the function is justified by the fact that, as \(1\) is the neutral element for multiplication in standard algebra (for all \(a\), \(1\times a = a\times 1 = a\)), \(e\) is the neutral element for multiplication in the max-plus algebra (for all \(a\), \(e \otimes a = a \otimes e = a\)).

xtilde_0 = zeros(18,1)
U = ones(2,5)

Trajectories \(\{\tilde{x}(k)\}_{k=0,1,\ldots,5}\) and \(\{y(k)\}_{k=0,1,\ldots,5}\) can be computed using functions trajectory and output_trajectory, respectively. Function trajectory is the same that we used in the previous section, but the input parameters are different here: trajectory(my_A, my_x0, k = my_k) computes the trajectory for an autonomous system, whereas trajectory(my_A, my_x0, B = my_B, u = my_u) computes the trajectory for a non-autonomous system. Functions trajectory and output_trajectory return matrices tr and out_tr, where each column corresponds to a state vector \(\tilde{x}(k)\) or output vector \(y(k)\) for different \(k\) values:

# compute state trajectory
tr = trajectory(A, xtilde_0, B = B, u = U)
# compute output trajectory
out_tr = output_trajectory(tr,C)

State trajectory:
ε     1.0     3.0     5.0     7.0     9.0
ε     3.0     5.0     7.0     9.0    11.0
ε     2.0     7.0    12.0    17.0    22.0
ε     7.0    12.0    17.0    22.0    27.0
ε    10.0    15.0    20.0    25.0    30.0
ε    22.0    27.0    32.0    37.0    42.0
ε       ε     1.0     3.0     5.0     7.0
ε       ε     3.0     5.0     7.0     9.0
ε       ε     2.0     7.0    12.0    17.0
ε ->    ε ->  7.0 -> 12.0 -> 17.0 -> 22.0
ε       ε    10.0    15.0    20.0    25.0
ε       ε    22.0    27.0    32.0    37.0
ε       ε       ε     1.0     3.0     5.0
ε       ε       ε     3.0     5.0     7.0
ε       ε       ε     2.0     7.0    12.0
ε       ε       ε     7.0    12.0    17.0
ε       ε       ε    10.0    15.0    20.0
ε       ε       ε    22.0    27.0    32.0

Output trajectory:
ε -> 22.0 -> 27.0 -> 32.0 -> 37.0 -> 42.0

Observe that the first six rows of tr correspond to the state trajectory \(\{x(k)\}_{k=0,1,\ldots,5}\).

Gantt charts

Gannt charts offer a way to visually represent the resource occupancy of places of a timed Petri net, which can be obtained from a state trajectory of the system. For this, we need to use the Gantt class and in the constructor we need the trajectory tr (already computed) as well as a list containing pairs of integers, which represent pairs of transitions. As an example for the above Petri net, if we are interested in the resource occupancy of the place between transitions \(x_1\) and \(x_2\) (i.e., the occupancy times of machine M1), we insert the pair (1, 2) into the list. For each pair of transitions, we will get one row in the chart showing the resource occupancy of the corresponding place as a function of time.

chart = Gantt(tr, [(1, 2), (3, 4), (5, 6)])

To obtain the figure, we can now call the draw method on this object. This takes one parameter which gives us the final time unit on the horizontal axis. More bars for the same resource stacked vertically indicate that more than one token is present in the place at a certain time.

chart.draw(40)

This figure can be saved as a pdf file, using the following command:

chart.save("figures/chart_name")

Figure saved as pdf: figures/chart_name

Precedence graphs

The library can analyze and draw precedence graphs of max-plus matrices. In the following code segment, we initialize a \(10 \times 10\) max-plus matrix and the precedence graph corresponding to this matrix.

A = Matrix([[  3,eps,  1,eps,eps,eps,eps,eps,eps,eps],
            [  3,eps,eps,eps,eps,eps,eps,eps,eps,eps],
            [eps,  4,eps,eps,eps,eps,eps,eps,eps,eps],
            [eps,  5,eps,  2,eps,eps,eps,eps,  2,  2],
            [eps,eps,  7,  1,eps,  3,  2,eps,eps,eps],
            [eps,  1,eps,eps,eps,eps,  2,eps,eps,eps],
            [eps,eps,eps,eps,eps,  2,eps,eps,eps,eps],
            [eps,eps,eps,eps,  1,eps,eps,eps,eps,eps],
            [eps,eps,eps,  4,eps,eps,eps,eps,eps,eps],
            [eps,eps,eps,eps,eps,eps,eps,eps,  1,eps]])

G = PrecedenceGraph(A)

The graph can be drawn using the draw method. It takes two optional arguments: critical, and spread. The first argument determines if we are drawing the critical graph or not (see later). The second one determines how far apart the nodes of the graph will be drawn. Here, we draw the graph (not the critical graph) with a spread of \(1\). Each maximal strongly connected component of the graph will have a different color, and gray nodes belong to no maximal strongly connected subgraph.

G.draw(False, 1)

We can see in the above figure that nodes of strongly connected components are grouped by color. Edges connecting them are of the same color.

We can for example obtain a list of maximal strongly connected components of the graph, a list of elementary cycles, and determine if the graph is strongly connected:

print("Maximal strongly connected components:", G.msc_components())
print("Cycles:", G.cycles())
print("Strongly connected:", G.is_strongly_connected())

Maximal strongly connected components: [[1], [4, 9]]
Cycles: [[1], [4], [1, 2, 3], [6, 7], [9, 4], [9, 10, 4]]
Strongly connected: False

The largest eigenvalue of the Matrix \(A\) is \(\lambda = 3\), which can be determined by this line of code:

A.eigenvalue()

3.0

The mean weight of a cycle is the sum of the weights of its arcs divided by the length of the cycle. The critical graph consists of cycles (and the corresponding nodes) which have a maximum mean weight (this equals to the largest eigenvalue \(\lambda\), in this case \(\lambda = 3\)). To draw the critical graph, we now use this line of code:

G.draw(True, 1)

We can see that each maximal strongly connected component is drawn in a different color. Nodes and arcs which are not part of the critical graph are colored gray and slightly transparent.

It is possible to save the matplotlib figures as a pdf file:

G.savefig("figures/precedence_graph")
G.savefig_critical("figures/critical_graph")

Saved precedence graph as pdf: figures/precedence_graph
Saved precedence graph as pdf: figures/critical_graph

References

• [1]: Hardouin, L., Cottenceau, B., Shang, Y., and Raisch, J., 2018. Control and state estimation for max-plus linear systems. Foundations and Trends® in Systems and Control 6, no. 1 (2018): 1-116.