%% Introduction to max-plus algebra
% Before proceeding, install the free Max-Plus Algebra Toolbox for Matlab from http://www.stanczyk.pro/mpa.
%% Scalar operations
% First create some data on which we are going to demonstrate the two max-plus algebra operations

a = 3;
b = 4;

% Max-plus addition is redefined according to a\oplus b = \max(a,b)

c = mp_add(a,b)

% Max-plus multiplication is redefined according to a\otimes b = a+b

d = mp_multi(a,b)

% There are also two special elements of the set. First, the max-plus zero, that is, a neutral element with respect to max-plus addition – oftentimes denoted \varepsilonin the literature

mp_zero

% and then the max-pluse one, that is, a neutral element with respec to the max-plus multiplication

mp_one

% which is sometimes denoted e in the literature, but we will not bother and use just 0.
% Now check that the two indeed work as neutral elements of the max-plus algebra with respect to the two operations. First, we compute a\oplus \varepsilon

mp_add(a,mp_zero)

% and then a\otimes 0

mp_multi(a,0)

%  Finally, the max-plus power

mp_power(a,3)

%
%% Operations with matrices
% First, define two matrices. They contain the \varepsilon = -\infty elements.

A = [2 3 mp_zero; 1 mp_zero 0; 2 -1 3]  
B = [mp_zero 5 -1; 3 mp_zero -2; mp_zero -4 7]

% We now demonstrate the same operations as we did for scalars above

C = mp_add(A,B)
D = mp_multi(A,B)
A2 = mp_power(A,2)

%% Computing the eigenvalue(s)

lambda = mp_mcm(A)

%
%% Demonstrating the result about eigenvalues of irreducible matrices
% Build the matrix A\oplus A^{\otimes^2 }

mp_add(A,A2)

% Since none of the off-diagonal elements is \varepsilon, the matrix is irreducible. This can also be checked by a dedicated function

mp_is_pgsc1(A)

% Then there exists integers c>0 and k\ge k_0, in this particular case

c = 1
k0 = 5

% such that for all k\ge k_0 we have A^{\otimes^{k+c} } =\lambda^{\otimes^c } \otimes A^{\otimes^k }
% Check it. We start with k = 4

mp_power(A,5) == c*lambda + mp_power(A,4)

% Obviously the condition is not satisfied. The reason is that k=4<k_0 =5. Now increase k

mp_power(A,6) == c*lambda + mp_power(A,5)
mp_power(A,7) == c*lambda + mp_power(A,6)

%
%% Solving Ax = B
% For the right-hand side

b = [1; 2; 3]

% the greatest subsolution is

x = mp_solve_Axb(A,b)

% We can check if what we got is really a solution

mp_multi(A,x) == b

% What we got is not the exact solution but rather the so-called greatest subsolution.
%
%% References
% De Schutter, Bart, and Ton van den Boom. "Model Predictive Control for Max-plus-Linear Discrete Event Systems." Automatica 37, no. 7 (July 1, 2001): 1049–56. https://doi.org/10.1016/S0005-1098(01)00054-1.
% ———. "Model Predictive Control for Max-plus-Linear Discrete-Event Systems: Extended Report & Addendum." Technical report. Control Systems Engineering. Delft, The Netherlands: Delft University of Technology, November 2000. https://pub.deschutter.info/abs/99_10a.html.
% De Schutter, Bart, Ton van den Boom, Jia Xu, and Samira S. Farahani. "Analysis and Control of Max-plus Linear Discrete-Event Systems: An Introduction." Discrete Event Dynamic Systems 30, no. 1 (March 1, 2020): 25–54. https://doi.org/10.1007/s10626-019-00294-w.