Tropical modeling of battery swapping and charging station

We propose and investigate a queueing model of a battery swapping and charging station (BSCS) for electric vehicles (EVs). A new approach to the analysis of the queueing model is developed, which combines the representation of the model as a stochastic dynamic system with the use of the methods and results of tropical algebra, which deals with the theory and applications of algebraic systems with idempotent operations. We describe the dynamics of the queueing model by a system of recurrence equations that involve random variables (RVs) to represent the interarrival time of incoming EVs. A performance measure for the model is defined as the mean operation cycle time of the station. Furthermore, the system of equations is represented in terms of the tropical algebra in vector form as an implicit linear state dynamic equation. The performance measure takes on the meaning of the mean growth rate of the state vector (the Lyapunov exponent) of the dynamic system. By applying a solution technique of vector equations in tropical algebra, the implicit equation is transformed into an explicit one with a state transition matrix with random entries. The evaluation of the Lyapunov exponent reduces to finding the limit of the expected value of norms of tropical matrix products. This limit is then obtained using results from the tropical spectral theory of deterministic and random matrices. With this approach, we derive a new exact formula for the mean cycle time of the BSCS, which is given in terms of the expected value of the RVs involved. We present the results of the Monte Carlo simulation of the BSCS's operation, which show a good agreement with the exact solution. The application of the obtained solution to evaluate the performance of one BSCS and to find the optimal distribution of battery packs between stations in a network of BSCSs is discussed.


Introduction
The usage of electric vehicles (EVs) has seen a great rise on a large scale in recent years [4,24].However, the adoption of electric vehicles is limited by problems such as the slow charging of battery packs (BPs) and the accelerated aging of BPs during fast charging (see, e.g., [13,6] for overviews of related problems and solution trends).Since batteries are the main source of power for EVs, ensuring energy supply is an important way to improve users' experience.At present, the power supply method for EVs is divided into two types: plug-in charging and battery swaps [15,13].Plug-in charging has disadvantages such as a long charging time, fast charging shortening the service life of the battery, and the parking lots required for charging taking up a larger space.In addition, if the daily load of residents and the peak of EVs' charging load are in the same time period, this will lead to a "peak-add-peak" state, which will affect the normal operation of the power grid.On the contrary, the battery swapping scenario addresses these problems well.Battery swapping for EVs can decrease user waiting time, reduce purchase cost, and improve batteries' useful life.Therefore, many companies have adopted the battery swapping scenario for EVs.However, there are still some challenges in the promotion of the battery swapping scenario for EVs, such as the operating cost of the battery swapping and charging station (BSCS) and the centralized battery charging load.
Modern research on the implementation of battery swapping offers a range of models to study various aspects of the BSCS's operations, including the battery logistics and transportation strategy as well as energy management and operation scheduling in EVs' battery swapping and charging systems and networks.A queueing model of BSCSs with Poisson arrival and constant service times was proposed by Choi and Lim in [5] to develop and analyze queuelength-dependent overload control policies.The queue distributions under different policies are derived using an embedded Markov chain, and system performance measures such as blocking probability and mean waiting time are examined by numerical examples.
In [1], the operation of a BSCS was represented using a finite-horizon Markov decision process model combined with a dynamic programming algorithm, which allows determining the number of BPs to recharge, discharge, and replace over time.In [36], an optimal scheduling problem was examined, which assigns the best BSCS to each EV, based on their current location and battery charge level.A scheduling strategy of the optimal transportation of BPs from a charging station to a swapping station was developed in [25].The strategy involves an optimization problem, which is solved using a genetic algorithm.This strategy was compared with two simple strategies by using the Monte Carlo simulation of battery swapping demand.
Given the cost and efficiency of operating battery replacement and charging stations, many scientists have proposed various strategies.Kang et al. [16] proposed a new strategy for the centralized charging of EVs within the framework of battery replacement scenarios.This strategy considers optimal charging priorities and charging locations and, ultimately, minimizes the total charging cost based on electricity spot prices.A battery-planning strategy based on a partitioned battery management method was developed by Yang et al. [34].In order to maximize profits, an optimization objective function was set, including the number of batteries in each segment.San et al. [29] obtained the optimal charging strategy for a single battery replacement and charging station in order to minimize the cost of charging it.The optimization model was transformed into a Markov decision process with constraints, and the optimal strategy was obtained using the Lagrangian method and dynamic programming.Liu et al. [26] demonstrated the method to optimize the charging and logistics of discharged and fully charged batteries to maximize the profits of battery replacement and charging stations.The problem of the optimal location of battery swapping stations (BSSs) was solved in [27] by using a fuzzy multi-criteria decision-making approach.
Wang et al. [32] addressed the problem of the online management of BSSs to minimize energy costs and ensure service quality.At the same time, the problem of designing optimal autonomous battery-swapping stations was studied to determine the optimal number of batteries and to achieve the ideal compromise between charging flexibility and battery cost distribution.Zhang et al. [37] confirmed that the peak demand for battery replacement services can be effectively reduced, and the use of optimization strategies can reduce the total cost by approximately 12%.The above research mainly considered the impact of factors such as the number of batteries, the battery charging costs, logistics planning, electricity prices during use, and the profits of battery replacement stations on the cost of battery replacement and charging stations to reduce operating costs.It suggested using different strategies to replace batteries and charging stations to varying degrees.However, the above view does not take into account the relationship between the load generated by large-scale centralized charging and the cost of replacing batteries and operating a charging station.
Huang et al. [14] found that the correspondence between the stochastic supply of wind and the need to charge electric vehicles can reduce the demand for conventional energy and carbon dioxide emissions.Xing et al. [33] proposed a data-based method for predicting electric vehicle charging demand based on online travel data.The predictive model provides recommendations for chargers and charging management.Battery charging load research mainly includes load prediction and the impact of reducing battery charging load on the power grid.The impact of reducing battery charging load is mainly focused on charging methods and increasing the overload capacity of the power grid.
However, there are few studies on the logistics and transportation of batteries at swapping and charging stations.Furthermore, current research does not propose specific transportation strategies to reduce the load generated by large-scale charging.The operating costs of battery replacement and charging stations and large-scale battery charging loads are the main areas for the implementation of battery replacement scenarios [31,35].
With the aim to contribute to the study of the BSCS's operations in this paper, the authors propose a framework based on tropical algebra to represent and analyze stochastic models of BSCSs.Tropical (idempotent) algebra deals with the theory and applications of algebraic systems with idempotent operations [2,18,7,12,8,21,3].A typical example of these systems is max-plus algebra, which is a semifield with addition defined as the operation of the maximum and multiplication as the arithmetic addition.One of the advantages of tropical algebra is that many problems, which are not linear in the ordinary sense, can turn into linear ones in the tropical algebra setting.The models and methods of tropical algebra find applications in various research domains such as location analysis, project scheduling, and decision-making.The application area includes stochastic dynamic systems, where tropical algebra serves as a useful tool to represent and analyze stochastic systems governed by tropical linear dynamic equations [12,11,21].
In this paper, we propose and investigate a queueing model of BSCSs for EVs.A new approach to the analysis of the queueing model is developed, which combines the representation of the model as a stochastic dynamic system with the application of the methods and results of tropical algebra.We describe the dynamics of the queueing model by a system of recurrence equations that involve random variables (RVs) to represent the interarrival time of incoming EVs.A performance measure for the model is defined as the mean operation cycle time of the station.Furthermore, the system of equations is represented in terms of the max-plus algebra in vector form as an implicit linear state dynamic equation.The performance measure takes on the meaning of the mean growth rate of the state vector (the Lyapunov exponent) of the dynamic system.By applying a solution technique of vector equations in tropical algebra, the implicit equation is transformed into an explicit one with a state transition matrix with random entries.The evaluation of the Lyapunov exponent reduces to finding the limit of the expected value of norms of tropical matrix products.This limit is then obtained using results from the tropical spectral theory of deterministic and random matrices.With this approach, we derive a new exact formula for the mean cycle time of the BSCS, which is given in terms of the expected value of the RVs involved.The numerical results of the Monte Carlo simulation of the BSCS's operation is demonstrated, which show a good agreement with the exact solution.We discuss the applications of the obtained solution to evaluate the performance of one BSCS and to find the optimal distribution of BPs between stations in a network of BSCSs.
As the obtained results show, the proposed approach based on the methods and results of tropical algebra presents a useful tool to model and analyze some classes of queueing models including the model of the BSCS under study.The approach can offer solutions, which are given in terms of the expected values of the RVs involved and independent of the particular probability distribution of the RVs.Such solutions may be of interest in the case when the details of the underlying probability distributions are difficult to determine and, thus, serve to supplement other modeling techniques with the need to fix a distribution.
The rest of the paper is organized as follows.In Section 2, we describe a queueing model of the BSCS that serves to both motivate and illustrate the study.Section 3 provides an overview of key definitions and the notation, and presents the preliminary results of tropical algebra, which are used in subsequent sections to examine the model under consideration.A stochastic dynamic model defined in the tropical algebra setting is described, and some related results are discussed in Section 4.An application of the tropical-algebra-based approach to the analysis of the BSCS model is demonstrated in Section 5, which includes a simple formula for calculating the mean cycle time of the BSCS and related simulation results.In Section 6, an example of the application of the obtained results to the optimal distribution of BPs between BSCSs is illustrated.Section 7 offers some concluding remarks.

Battery Swapping and Charging Station Model
We considered a BSCS that serves incoming requests of EVs to swap a depleted (discharged) BP to a fully charged one.Each EV is assumed equipped with one BP, and all BPs are considered of the same type (identical).The BSCS consists of a battery swapping and battery charging/storage areas.
The station has a set of identical BPs located in the storage area where they are charging and, then, waiting for use in swapping.All BPs can be charged simultaneously, and the charging of every BPs takes the same time.The swapping operations are performed one at a time and require equal time for all EVs.
The EVs arrive at the BSCS at random with time intervals distributed according to some probability law.Upon arrival, an EV waits until the following conditions hold: (i) a fully charged BP is available, and (ii) the swapping of the BP for the previous EV is completed, or the swapping procedure is immediately started if a fully charged BP and the swapping unit are both available.
A graphical representation of a BSCS as a queueing model is given in Figure 1.The model consists of (i) a single-server queue, which represents an arrival source of EVs, (ii) a single-server fork-join queue, which represents the swapping of batteries, and (iii) a multi-server queue, which represents the charging of the batteries.All queues have infinite buffers.At the initial time, the first queue is assumed to have an infinite number of jobs (EVs), the second queue has no job (the EV and BP ready for swapping), and the third queue has m jobs (BPs ready for charging).

System of Recurrence Equations
Suppose that the station has m BPs intended for swapping, and define the following state variables.For k = 1, 2, . .., let x(k) be the arrival epoch of the kth incoming EV, y(k) be the completion time of the battery swapping for the kth EV, and z(k) be the time when a fully charged BP is available for the kth EV.We also assumed that x(k) = y(k) = z(k) = 0 for all k ≤ 0. We now describe the evolution of the system as a set of recurrence equations.We denote the time interval between the (k − 1)the and kth arrival epochs by α k and assume {α k | k = 1, 2, . ..} to be a sequence of independent and identically distributed positive (nonnegative) RVs with a finite expected value Eα 1 = a ≥ 0 and variance Dα 1 .With this notation, we can represent the kth arrival epoch as Furthermore, we denote by b > 0 the swapping time of one BP.Observing that the kth swapping operation starts as soon as the following events occur: (i) the kth EV arrives, (ii) the (k − 1)the swapping operation completes, and (iii) a fully charged BP becomes available for the kth time, we write the equation: Finally, we assumed that all m BPs available in the station are discharged at the initial time epoch k = 0.With the charging time of one BP denoted by c > 0, we have We now substitute z(k) from the last equation into the second and, then, combine the first and second equations into the dynamic system in two state variables: (1)

Performance Measure
We define the operation cycle of the BSCS model described by (1) as the interval between successive completions of swapping operations.Furthermore, we consider the mean (average) cycle time over the first k cycles: We turn to the limit when k tends to ∞ and assume that this limit exists (deterministically or with probability one) to write The constant λ is referred to as the mean cycle time and may serve as a useful characteristic of the system.Specifically, for a large time horizon T , the ratio T /λ differs little from the mean number of battery swaps in the time interval from 0 to T .Since this ratio also shows the mean number of batteries swapped, it can be used to estimate other characteristics such as the mean total energy consumption for battery charging (which is considered proportional to the mean total charging time estimated as cT /λ, where c denotes the energy consumption per BP) or the mean total revenue received from customers (proportional to T /λ).
The evaluation of the mean cycle time directly from recurrent Equation (1), which involves the operation of the maximum, may be a rather difficult problem even if the equations are in a simple form like those presented above.At the same time, this form of the recurrence equations offers a potential for the use of models and methods of tropical algebra, which allows one to handle such equations in a unified analytical framework.In the subsequent sections, we show how to represent the equations in terms of tropical algebra in vector form and, then, use this representation to evaluate the mean cycle time analytically.

Elements of Tropical Algebra
We begin with preliminary definitions and results of tropical algebra, which are used for the representation and analysis of the dynamic model in what follows.Tropical (idempotent) algebra deals with the theory and applications of algebraic systems with idempotent operations, which are studied in many works, including the monographs [2,18,7,12,8,21,3].

Idempotent Semifield
Let X be a set that is closed under associative and commutative binary operations: addition ⊕ and multiplication ⊗, and includes their neutral elements: zero 0 and identity 1. Addition is idempotent: x ⊕ x = x for all x ∈ X. Multiplication distributes over addition and is invertible: for each x = 0, there is an inverse x −1 such that x ⊗ x −1 = 1 (hereafter, the multiplication sign ⊗ is omitted to save writing).
The power notation with integer exponents specifies iterated products: x p = x p−1 x, x −p = (x p ) −1 , 0 p = 1, and x 0 = 1 for x = 0 and integer p > 0. The powers with rational exponents are also assumed well-defined.The binomial identity takes the form of the equality (x⊕y) q = x q ⊕y q , which is valid for any rational q ≥ 0.
The set X is assumed to be totally ordered by an order relation consistent with that induced by idempotent addition by the rule: x ≤ y if and only if x ⊕ y = y.With respect to this order, addition and multiplication are monotone in each argument: if the inequality x ≤ y holds, then x ⊕ z ≤ y ⊕ z and xz ≤ yz for any z.For nonzero x and y such that x ≤ y and rational q, the inequality x q ≤ y q holds if q ≥ 0 and x q ≥ y q if q < 0. Furthermore, the inequalities x ≤ x ⊕ y and y ≤ x ⊕ y are valid for all x and y.Finally, the inequality x ⊕ y ≤ z is equivalent to the system of inequalities x ≤ z and y ≤ z.
A typical example of the system is the real semifield (R ∪ {−∞}, max, +, −∞, 0), also known as max-plus algebra.In max-plus algebra, the operations are defined as ⊕ = max and ⊗ = + and the neutral elements as 0 = −∞ and 1 = 0.For any x ∈ R, the multiplicative inverse x −1 is equal to the opposite number −x in the standard arithmetics.The power x y coincides with the usual arithmetic product x × y.The order relation ≤ corresponds to the natural linear order on R.

Algebra of Matrices and Vectors
The scalar operations ⊕ and ⊗ are extended to vectors and matrices over X in the usual way.A matrix with all entries equal to 0 is the zero matrix denoted 0. For any matrices A = (a ij ), B = (b ij ) and C = (c ij ) of appropriate sizes, and scalar x, matrix addition, matrix multiplication, and scalar multiplication are defined by componentwise formulas: For any nonzero (m × n)-matrix A = (a ij ), its multiplicative conjugate is the (n × m)-matrix The monotonicity of scalar addition and multiplication, as well as other properties that involve the order relations are extended to the operations on matrices, where the inequalities are understood componentwise.
A square matrix is diagonal if all its off-diagonal entries are equal to 0 and triangular if its entries either above or below the diagonal are equal to 0. A triangular matrix with all diagonal entries equal to 0 is called strictly triangular.The block diagonal and (strictly) block triangular matrices are introduced in a similar way.
A diagonal matrix that has all diagonal entries equal to 1 is the identity matrix denoted by I.The power notation for matrices is defined in the sense of tropical algebra as follows: A 0 = I, A p = AA p−1 for any square matrix A and integer p > 0.
The trace of a square matrix A = (a ij ) of order n is given by A tropical analogue of the matrix determinant is defined as If the condition Tr(A) ≤ 1 holds, the Kleene star matrix is calculated as A matrix that consists of one column (row) is a column (row) vector.All vectors are assumed column vectors unless transposed.A vector with all entries equal to 0 is the zero vector denoted 0. Any vector that has no zero entries is called regular.The vector that has all entries equal to 1 is denoted by 1 = (1, . . ., 1) T .In max-plus algebra, where 1 = 0, the vector 1 has all entries equal to arithmetic zero (the usual zero vector).
For any matrix A = (a ij ) and vector x = (x i ), tropical norms are given by which coincide in max-plus algebra with the maximum entries of A and x.
For any conforming matrices A, B and C, and scalar x, the following relations hold: A scalar λ is an eigenvalue of an (n × n)-matrix A if there exists an n-vector x = 0 such that Ax = λx.The spectral radius of A is the maximum eigenvalue, which is given by Note that, if the spectral radius is defined in the framework of max-plus algebra, it can be represented using ordinary arithmetic operations as the maximum of the mean (average) cyclic sums of entries in A in the form: If a matrix A has no entries equal to 0, then for all integers k ≥ 0, the following inequality holds (see, e.g., [23]): The next theorem is a consequence of the results obtained in [30,28] (see also [18,23]).
theorem 1.For any (n × n)-matrix A, there exist the limits:

Vector Equation and Matrix Inequality
To conclude the overview of the preliminary results, we give a solution for a vector equation and derive inequalities for products of square matrices to be used in what follows.Suppose that given an (n × n)-matrix A and n-vector b, the problem is to find regular n-vectors x that satisfy the equation: The following lemma offers a solution of the equation in a special case as a consequence of the results obtained in [19,21].We now turn to evaluating the lower and upper bounds for the norm of a product of matrices in block triangular form.Let A(i) for all i = 1, . . ., k be conforming block triangular matrices given by the sum of block diagonal and strictly triangular matrices as follows: Consider the product of the matrices A(i) over all i = 1, . . ., k, and denote it by To simplify further formulas, we introduce the notation: where the empty products are thought of as equal to the identity matrix I.
The next statement offers lower and upper bounds on the norm A k .
proposition 3. Let A(i) for all i = 1, . . ., k be matrices defined as (6).Then, the following double-inequality holds: Proof.To obtain a lower bound, we use the inequality A(i) = D(i) ⊕ T (i) ≥ D(i), which holds for all i = 1, . . ., k.By combining these inequalities, we have Furthermore, we consider the upper bound.The application of the distributivity of multiplication over addition yields We consider the product under summation and apply the properties of the norm to write for all i, we obtain the upper bound: We note that the matrix D k is block diagonal, and hence, D k = D 1k ⊕ D 2k .It remains to combine both the lower and upper bounds, which yields (7).

Stochastic Dynamic Systems
In this section, we examine stochastic dynamic systems in the tropical algebra setting, which are used for the description of the evolution of the queueing system under study in the next section.The main purpose of this section is to evaluate the Lyapunov exponent for a dynamic system with a state transition matrix of special form.For further details on the application of tropical algebra to stochastic dynamic systems, one can consult [12,11].
We consider a dynamic model that is governed by the state equation represented for all k = 1, 2, . . . in terms of max-plus algebra in the form: where x(k) denotes a state n-vector and A(k) a state transition (n × n)-matrix given by Each entry a ij (k) of the matrix A(k) may be an RV or a constant.The corresponding random entries in the matrices A(k) for k = 1, 2, . . .are assumed independent and identically distributed (i.i.d.) with finite expectation and variance.Note that the random entries in one matrix A(k) need not be independent.
We define the matrix product: With this notation, the state dynamic equation at (8) can be reduced to The Lyapunov exponent indicates the mean growth rate of the state vector, and it is defined as the limit: We note that, in the context of max-plus algebra, the last definition is represented in the conventional form: Furthermore, with x(0) = 1 where 1 = (0, . . ., 0) T , we have and then, rewrite the above limit as The next result from [20,21] (see also [12,11]), which is a consequence of Kingman's subadditive ergodic theorem [17], gives conditions for the above limit to exist with probability one (w.p. 1) and evaluates λ as the limit of the expected values of A k 1/k .theorem 4. Let {A(k)| k = 1, 2, . ..} be a stationary sequence of random matrices such that Then, there exists a finite number λ such that Since the matrices A(k) are assumed to be i.i.d., the sequence of these matrices is stationary.Moreover, since the random entries in A(k) have finite expected values, the condition E A 1 < ∞ holds.The condition ρ(E[A 1 ]) > −∞ actually means that the sequence of matrices A k does not degenerate into a zero matrix 0 (which has all entries equal to −∞ in max-plus algebra), and it is assumed satisfied.
It follows from Theorem 4 that, for the dynamic systems under consideration, the Lyapunov exponent exists and can be found as the limit of the expected values E A k 1/k as k tends to ∞.The evaluation of the limit and of the expectations E A k themselves can be a difficult problem.However, it is not difficult to solve the problem for matrices A(k) that have a particular form or structure [20,21].Specifically, if the matrices A(k) are triangular, then the Lyapunov exponent is calculated as In the context of max-plus algebra, the above formula turns into the maximum of the expected values of the diagonal entries in A 1 = A(1) given by λ = max 1≤i≤n E[a ii (1)].
We note that the same result is valid for the diagonal matrices A(k) as well.Moreover, this result can be readily extended to the system (8) with block diagonal matrices.lemma 5. Let A(k) for k = 1, 2, . . .be block diagonal matrices of the form: Consider matrices , and suppose that D rk 1/k → µ r w. p. 1 as k → ∞ for all r = 1, . . ., s.Then, the Lyapunov exponent of the system (8) is given by λ = s r=1 µ r .

we apply the binomial identity to write the equality
. It remains to let k go to ∞ on both sides of the equality, which yields the desired result.
The extension of this result to block triangular matrices seems to be not so easy.Below, we evaluate the Lyapunov exponent for block triangular matrices of special form.Consider a dynamic system with state transition matrices A(k) of block triangular form defined as (6) in the framework of max-plus algebra.We suppose that the upper diagonal block reduces to an RV α k and the lower is a constant nonrandom matrix D to write We assume that α k for k = 1, 2, . . .are i.i.d.RVs that have a finite expected value and variance, and the matrix D has no zero entries.The RVs T (k) are also assumed i.i.d. with nonnegative expectation and finite variance.lemma 6.Let Eα k = µ 1 ≥ 0 be the expected value of α k and ρ(D) = µ 2 > 0 be the spectral radius of D. Then, the Lyapunov exponent of the system is given by λ Proof.To verify the statement, we show that (7), which yields First, we examine the right inequality.We apply (4) we can write the inequalities: With these inequalities, we expand the right inequality at (9) as follows: Next, we rewrite the last inequality in terms of ordinary operations and take expectations.Taking into account that E DD − = DD − , we obtain We note that DD − is bounded, and hence, DD − /k → 0 as k → ∞.Furthermore, T (j) for j = 1, 2, . . .are i.i.d.RVs with finite expectation and variance.As k goes to ∞, the expected value of the maximum of these RVs grows as O(k 1/2 ) [9, 10], and therefore, Consider the last term on the right-hand side of (10), and suppose that Eα 1 = µ 1 ≤ µ 2 .We represent this term as We observe that α i − µ 2 are i.i.d.RVs with the expectation E(α i − µ 2 ) ≤ 0 and finite variance.Since the expected value of the maximum of the cumulative sums of these variables grows as O(k 1/2 ) as k tends to ∞ (see, e.g., [22]), we have Using similar arguments, we can verify that if µ 1 ≥ µ 2 , then As a result, we conclude that Consider the left inequality at (9).As k tends to ∞, we have in terms of the usual operations, we see that Therefore, the left inequality leads to the inequality: Since the opposite inequality holds, we arrive at the conclusion that which completes the proof.
It is not difficult to see that this result remains valid if the matrix D may have zero entries, but some of its power D p is a matrix without zero entries.Indeed, in this case, we can consider a dynamic system: x where we use the notation: The matrix A ′ (k) has a block triangular form with D p as its lower diagonal block.

Application to Battery Swapping and Charging Station Model
We are now in a position to apply the previous results to represent the BSCS queueing model in terms of max-plus algebra and evaluate the mean operation cycle time for the model.We start with scalar recurrence Equations (1), which describe the dynamics of the model, and represent these equations in terms of max-plus algebra.Next, we rewrite the scalar equations in vector form as an implicit vector equation.This equation is solved to obtain an explicit state dynamic equation with a state transition matrix having random entries.Finally, by applying results from the previous sections, we derive an exact formula for calculating the mean cycle time under evaluation.

Tropical Representation of the Model
Let us rewrite the equations in (1) in terms of max-plus algebra.After replacing the operation max by the addition ⊕ and + by the multiplication ⊗ (the sign ⊗ is eliminated from the subsequent expressions), the equations become linear in the tropical sense and take the form: To represent the dynamic system in vector form, we introduce the following vector and matrices (where we use the notation 0 = −∞ and 1 = 0): . . .
With this notation, the system is written as an implicit equation in v(k) in the form: We solve this equation for v(k) by using Lemma 2. First, we note that tr B = 0. Furthermore, we see that B 2 = 0, and hence, B i = 0 for all i ≥ 2. As a result, we have Tr(B) = 0 and calculate The application of Lemma 2 leads to the explicit state dynamic equation: with the state transition matrix:

Tropical Representation of Performance Measure
We now exploit the dynamic model derived above to evaluate the mean cycle time λ given by (2).First, we see from scalar Equations (1) that the following inequalities are valid: As a result, we obtain y(k) = max(x(k), y(k), . . ., y(k − m + 1)).
Since the right-hand side of the above equality coincides with the max-plus algebra norm v(k) , we conclude that Therefore, the mean cycle time (2) can be represented in terms of max-plus algebra as By Theorem 4, we can find the mean cycle time by evaluating the limit of the expected values as follows:

Evaluation of Mean Cycle Time
To evaluate the mean cycle time of the system, we apply Lemma 6.Consider the state transition matrix A(k) at (11), and note that it has the block triangular form: where the matrix blocks are given by As is easy to see, the state transition matrix A(k) has the same form as in Lemma 6 and satisfies the assumptions of this lemma.Moreover, it is not difficult to verify by direct computation that the matrix D m−1 has no zero entries.
It follows from Lemma 6 that the Lyapunov exponent (the mean cycle time) of the system is given by λ The evaluation of the spectral radius ρ(D) by using (3) yields As a result, the mean cycle time is represented in terms of max-plus algebra as After rewriting in terms of the conventional algebra, we have Finally, we note that the obtained formula takes into account the expected value a = Eα 1 of the random interarrival time of incoming EVs and does not require a complete description of the underlying probability distribution.To illustrate this result by numerical experiments, we used Monte Carlo simulation to estimate the mean cycle time of the system over a long time horizon.
The simulation model consists of the recurrence relations in (1), which formally describe the evolution of the system.The simulation experiment involves the calculation of the state variables x(k) and y(k) for k = 1, . . ., K, where K is sufficiently large.The calculation of x(k) includes sampling from a given probability distribution to fix an interarrival time for each k.The estimates λ(k) = y(k)/k of the mean cycle time are evaluated for successive k to observe the convergence of the estimates to the value of λ specified by (12).
In Figures 2-5, we demonstrate the results of estimating the mean cycle time λ for a BSCS model with the following parameters fixed: the swapping and charging time of one BP were set to b = 5 and c = 100 and the number of BPs available at the station to m = 4.As an interarrival time distribution for EVs, the exponential and uniform distributions with means a = 25 and a = 30 were considered.The number of EVs was set to K = 200.
Figures 2 and 3 show the results of evaluating the estimators λ(k) = y(k)/k for k from 1 to 200 step 5 by the simulation of the system with the exponential interarrival time with mean a = 25 and a = 30, respectively.The simulation results for the same BSCS with the interarrival time uniformly distributed over [5,45] (with mean a = 25) and [10,50] (with mean a = 30) are given in Figures 4 and 5  be explained by the high arrival rate of EVs, which makes, in this case, the BP recharging and swapping process a dominated time factor that diminishes the influence of random fluctuations in the interarrival time.

Example of Application Problem
In this section, we offer an example of the application of the obtained results to solve real-world problems.Consider a network that consists of N BSCSs.For each station i = 1, . . ., N , let a i be the mean interarrival time of EVs.We denote the swapping time and the charging time of one BP by b i and c i , respectively.
Assume that the BSCS is equipped with m i BPs intended for swapping, and examine the mean cycle time for station i, which is given by The mean swapping rate at the station is evaluated as 1/λ i , whereas the mean number of BPs swapped for a large time horizon T is T /λ i .
Let us suppose that one swapping at station i generates an income r i .Then, the mean total income during time T is equal to .
We represent the mean total income as a function of the number m of BPs at the station in the form It follows from the representation that the function R i (m) increases until m becomes greater than a threshold value (b i + c i )/ max(a i , b i ), and remains unchanged with further increase of m.The maximum mean total income and corresponding optimal number of BPs are defined as , where [x] denotes the integer part of x.
Suppose there are M BPs, which we need to distribute between the BSCSs in the network so as to minimize (maximize) an appropriate optimality criterion.If the purpose is to maximize the mean total income generated by the network, the problem is formulated to find the number m i of BPs for each station i to attain the maximum: As a reasonable approximate solution technique, we can define the optimal numbers m i to be proportional to w i = r i /(b i + c i ).With this technique, the number m i is first found for each i = 1, . . ., N as the nearest positive integer: Furthermore, we check whether the numbers m i are outside their threshold values or not.If, for each i, the inequality m i ≤ (b i + c i )/ max(a i , b i ) holds, then the obtained numbers m i are taken as a solution to the problem.
Suppose that m i > (b i + c i )/ max(a i , b i ) for some i.In this case, we decrement m i by one and increment some m j such that j = arg max We continue to redistribute BPs between stations until all stations have the number of BPs within their threshold values.

Conclusions
In this paper, we propose a new approach to the analysis of BSCSs' operation, which combines queueing modeling with the application of the methods and results of tropical algebra.We started with the development of a queueing model in the form of a system of recurrence equations that determine the dynamics of a BSCS.We introduced a related performance measure in the form of the mean operation cycle time.Then, the model was represented in terms of max-plus algebra as a linear vector dynamic system with a random state transition matrix, whereas the performance measure became the Lyapunov exponent of the system.We applied the methods and techniques of tropical algebra together with the results on the convergence of the expected value of the maximums of random variables to find the Lyapunov exponent as a limit of the expected value of the matrix norms.After the calculation of the Lyapunov exponent, we arrived at an explicit expression in terms of the expected values of the random variables and constants involved.We showed how this expression can be used to evaluate and optimize the performance of the BSCSs' operation.
We believe that the described research demonstrates the strong potential of the proposed approach to investigate various dynamic models that can be represented as stochastic linear dynamic systems in the tropical algebra setting.The results obtained indicated the ability of the approach to supplement and complement existing techniques of the modeling and optimization of BSCSs' operation.The approach offers a potential to provide explicit results that are given in terms of the expected values of the the RVs involved in the dynamic model and do not require a specific probability distribution.At the same time, the approach can be applied only to dynamic models that are linear in the tropical algebra sense such as queueing models, whose dynamics can be described in terms of the operations of the maximum and addition.This constitutes one of the limitations of the obtained solution, which may make it difficult to extend this approach to other classes of queueing models.Another limitation is that the approach focuses on evaluating the Lyapunov exponents of dynamic systems and can be hardly extended to other performance measures of interest.
Possible directions of future research may concern further investigation of the obtained solution, including the sensitivity analysis of the model.An extension of the BSCS model to incorporate more-complicated operation patterns and accommodate additional constraints are of particular interest.As an example, one can consider a station where the number of simultaneously charged BPs is limited or the battery charging time is random.The formulation of new meaningful optimization problems to improve BSCSs' performance and the development of efficient solutions constitute another promising line of investigation.

Figure 1 :
Figure 1: Queueing model of battery swapping and charging station.

25 Figure 2 : 30 Figure 3 :
Figure 2: Simulation results for the exponential distribution with expected value of 25.