Koopman Mode Analysis of agent-based models of logistics processes

Modern logistics processes and systems can feature extremely complicated dynamics. Agent Based Modeling is emerging as a powerful modeling tool for design, analysis and control of such logistics systems. However, the complexity of the model itself can be overwhelming and mathematical meta-modeling tools are needed that aggregate information and enable fast and accurate decision making and control system design. Here we present Koopman Mode Analysis (KMA) as such a tool. KMA uncovers exponentially growing, decaying or oscillating collective patterns in dynamical data. We apply the methodology to two problems, both of which exhibit a bifurcation in dynamical behavior, but feature very different dynamics: Medical Treatment Facility (MTF) logistics and ship fueling (SF) logistics. The MTF problem features a transition between efficient operation at low casualty rates and inefficient operation beyond a critical casualty rate, while the SF problem features a transition between short mission life at low initial fuel levels and sustained mission beyond a critical initial fuel level. Both bifurcations are detected by analyzing the spectrum of the associated Koopman operator. Mathematical analysis is provided justifying the use of the Dynamic Mode Decomposition algorithm in punctuated linear decay dynamics that is featured in the SF problem.


Introduction
An agent-based model (ABM) is a computational technique in which behavior of individual agents is encoded by simple rules, and the outcomes are observed at the scale of the system [1]. Simple rules of behavior for individuals can lead to complex, system-scale emergent phenomena [2,3]. ABMs are relevant to modern modeling problems because their essential feature of interacting micro-scale agents leading to macro-scale dynamics resembles many of the decentralized, highly-interacting social [4][5][6][7][8][9] and economic [10,11] systems of today.
Beyond demonstrating the existence of emergent phenomena for their own sake, the prime importance of ABM results is their analysis for the purpose of decision making. Agent-based models have been analyzed with applications to both governmental and business decision making interests. Das and Hanaoka [12] examined an ABM of humanitarian logistics applied to the allocation of relief supplies by disparate stakeholders. Démare et al [13]  (DMD) algorithm used to compute Koopman eigenvalues and modes, then introduce the MTF and SF logistics simulations. In the Results section we describe the results of the MTF and SF logistics simulations and the understanding of their dynamics made possible by KMA. We also provide a theoretical description of KMA applied to linear decay dynamics with random increases and show that KMA theory correctly predicts a characteristic result observed in the analysis of ABM results. The Discussion and Conclusions section discusses the significant findings of this work and presents our conclusions about these systems and the utility of KMA in describing them.

Koopman Mode Decomposition (KMD)
Koopman Mode Analysis (KMA) decomposes a complex signal, fg t 2 R n : t 2 I � Rg, into a linear combination of spatial structures that have simple temporal dynamics plus a "noise" term as It does this by computing a data-driven spectral analysis of an underlying family of linear operators which induce the evolution of the signal. In (1), λ j = a j + iω j is an eigenvalue of the underlying family of operators, a complex frequency determining the temporal signature; v j is the mode which determines the spatial structure of the associated temporal signature; the overlines denote complex conjugation; and n t 2 R n is the "noise" term due to the continuous spectrum of the family of operators. We remark that the noise term does not necessarily need to have Gaussian statistics as is assumed in many other contexts.
The family of operators is induced as follows. We assume that fg t 2 R n : t 2 Ig is a stochastic process. There is an underlying state space Z, a flow S t : Z ! Z (t 2 I), and a vector-valued function g : Z ! R n such that The composition operation on the right-hand side induces a family of linear, infinitedimensional operators (called the family of Koopman operators) U t : F ! F , where F ¼ fg : Z ! R n g is a vector space of functions, satisfying and therefore the original signal can be written in terms of the family of Koopman operators as Since {U t } is a family of linear operators, we can analyze its spectrum. Assuming its a spectral family, the family can be decomposed as U t g ¼ X j e l j t v j ðgÞ þ e l j t v j ðgÞ þ n t ðgÞ: Using this with (4), gives (1) (where we have suppressed the dependence of the v j 's and n t on g).
In KMA, the eigenvalues, e l j t , and the modes, v j , from (4) are determined via the spectral analysis of a data-driven approximation of the family of Koopman operators U t .
In terms of a logistics problem such as the Medical Treatment Facility (MTF) problem, the components of g t could represent the number of patients at each of the n medical facilities at time t. The state space and flow map (Z and S t , respectively) would correspond to the state space for the underlying ABM model and its simulation, respectively. A Koopman mode v j ¼ ðv ð1Þ j ; . . . ; v ðnÞ j Þ T would determine which of medical facilities had its level of patients vary with a temporal signature partially characterized by the Koopman eigenvalue λ j = a j + iω j . For example, if v ðiÞ j was non-zero, then this would indicate that the number of patients at facility i is partially characterized by growth (a j > 0) or decay (a j < 0) with rate e a j and oscillation with frequency ω j . For any fixed medical treatment facility i, its full temporal behavior would be given by where g ðiÞ t , v ðiÞ j , and n ðiÞ t are the i th components of g t , v j and n t , respectively. The most popular computational algorithm for computing KMA's eigenvalues and modes is the Dynamic Mode Decomposition (DMD) algorithm which has been solidified as an important tool in the data-driven analysis of complex dynamical systems. First introduced into the Fluids community by P. Schmid in 2008 [20,21]-although without reference to Koopman operator theory-there has been much subsequent research effort in improving the basic algorithm's, and its variants', correctness, stability, computational efficiency, and theoretical underpinnings. The book [22] gives a good introduction to the various methods and algorithms. In this paper, we use a variant of DMD recently introduced in Drmać et al. ([23], Algorithm 2) which includes an optional refinement procedure for the computed modes and eigenvalues as well as an explicit error term specifying the accuracy of the computed spectrum. We discuss the basic DMD algorithm here and refer the reader to [23] for details on the refinement procedure.
The DMD algorithm requires a finite sequence of snapshots of the signal, (f 0 , . . ., f m ), where f k = g kΔt , as an input and returns the matrix V = (v 0 , . . ., v m−1 ) and the vector Γ = (γ 0 , . . ., γ m−1 ). The i-th column of V and the i-th element of Γ are approximations of some v j from (1) and the associated eigenvalue e l j Dt , respectively. Note that f k = U Δt f k−1 . The DMD algorithm first looks for a matrix A : R n ! R n that approximates U Δt on the Krylov subspace K ¼ spanff 0 ; . . . ; f mÀ 1 g � F . That is, we are looking for a matrix representation A of P K U Dt : K ! K, where P K : F ! K is the orthogonal projection of functions onto K. The DMD algorithm then uses this matrix to approximate the Koopman modes and eigenvalues, v j and e l j Dt , respectively.
Algorithm 1: DMD Remark. Note, we are indexing from 0 in Algorithm 1, K k is an orthogonal basis for the kdimensional Krylov subspace K, and A is the matrix representation of P K U Dt : K ! K in the SVD basis K k . The Koopman eigenfrequency λ j can be computed from γ j via l j ¼ 1 Dt log ðg j Þ, where we take the principal value of the complex logarithm.

KMA and logistics problems
In terms of the logistics problems, f k 's are quantities of interest during the evolution of the agent-based model of the logistics problem-e.g., the amount of fuel at each of the sites. At a very coarse-level the use of KMA/DMD for these problems is a recurrent 3-step procedure as time evolves.
1. Measure quantities of interest, f k , at time t k (either from the simulation of an underlying ABM model or from a real system).
2. Input finite windows of snapshots (f k−m , . . ., f k ) into the DMD algorithm and return Koopman modes and eigenvalues v j and λ j for that time window.
3. Use the eigenvalues and modes to interpret the overall health of the system.

Update simulation time and goto step 1:
Step 3 is the most involved step and is problem dependent. However, our goal is to use the Koopman spectrum to interpret the overall health of the systems under observation. In the results below, there will be a change in the spectrum as a simulation parameter is changedthe casualty source rate in the MTF problem, the maximum fuel capacity of a site in the ship refueling problem. As the eigenvalue distribution changes with changes in the parameter, the system will undergo a "bifurcation", analogous to how the stable spiral in a Hopf bifurcation changes into an unstable spiral and a stable limit cycle appears as the eigenvalues of the fixed point cross the imaginary axis.
In the logistics problems, there do not exist obvious events to use in determining the bifurcation, such as there is with the Hopf bifurcation. There is a possibility of developing such theory for stochastic systems using an operator-theoretic version of spectral stability theory analogous to the one developed in [24]. The basic ingredients necessary for such theory have been developed for stochastic systems in [25], but are not available as full stability and bifurcation theory yet. Thus, here we take a more heuristic approach. In the results below, we will see the eigenvalues clustering around 0 (or 1) and then moving toward 1 (or 0) as the parameters are changed. While the distributions of the eigenvalues at the extremes of the parameter values is easy to see visually-and we can easily say the system behaves qualitatively differently for these parameter values-defining a unique bifurcation point is problematic. One way to automatically do this is as follows. Define a circle of radius r, centered at 1 in the complex plane, and define a threshold value 0 < h � 1. For each parameter value, θ, compute the percentage of Koopman eigenvalues lying inside the circle around 1 One definition of the bifurcation point is to define it as the parameter value, θ � , for which This definition, however, is not particularly useful for a stochastic system since θ � is not unique. A better way is to define two thresholds 0 < h g < h r � 1. We then say the system is healthy ("in the green") if and unhealthy ("in the red") if This gives operators an automatic method of quickly interpreting the health of the system. Of course, the usage of "healthy" and "unhealthy" above is dependent upon the system actually exhibiting undesirable behavior if the eigenvalues cluster around 1; for other systems this may be desirable behavior. Additionally, centering the r-circle around a different eigenvalue may make more sense for a particular system. Finally, the values of the parameters r, h g , h r must be determined for each system under study according to the desired performance characteristics. These are topics that will be addressed in future research.

Medical treatment facility logistics simulation overview
The medical treatment facility (MTF) model consists of a network of 15 medical treatment facilities and five sources representing combat zones that produce casualties. Each time step in the simulation represents one day. Each source generates the same number N c of casualties per day, where N c = 30, 50, . . ., 330 depending on the simulation. For each MTF and source pair, there is a given probability that a casualty would be directed from that source to that MTF. The distance between each MTF and source is fixed for each pair, and for each kilometer of distance between the pair a casualty has a 1.5% chance of becoming Dead of Wounds (DOW), where each kilometer traveled is a Bernoulli trial. Each MTF has a patient capacity, where for each simulation the sum of the MTF capacities equals 800 casualties but the capacity of each MTF is randomly determined, subject to the constraint that no MTF capacity could be less than 5 casualties. Each MTF also has a fixed probability that, for each time step (day) in the simulation, a patient would become DOW or return to duty (RTD), and also an output rate, representing the number of patients (up to the occupancy of the MTF) to be transferred from the MTF to a hospital and exiting the simulation. At the end of each time step, all patients in an MTF in excess of the MTF's capacity become DOW. Fig 1 shows the connectivity of the MTF network, where each arrow indicates a connection between each source to each MTF. Table 1 shows the MTF to source distance for each pair and Table 2 shows the probability that a casualty from that source will be directed to a given MTF. Table 3 shows the DOW and RTD probabilities and output rates for each MTF.

Ship fueling logistics simulation overview
The second model is an abstract representation of the logistical supply of fuel to fixed-site facilities by transport ships. The facilities are represented by a static network of 36 mission sites on a 6-by-6 hexagonal grid (see Fig 2), and the transport ships are represented by six dynamic agents referred to as assets. There is an additional site, referred to as the depot, which is distinct from the mission sites. It is connected to two of the mission sites and serves as an infinite source of fuel units.
Fuel is necessary for sites and assets to function. Each site has a single fuel level, and when the fuel level of any site reaches zero the simulation ends, with the number of time steps completed called the "mission life" for that simulation. All sites begin the simulation with an equal number of fuel units, which in this work was varied from 25 to 200 units in steps of 25. Assets have two fuel levels: 1) an internal level of the asset itself, which has a starting and maximum value of 150 fuel units and if equal to zero results in the loss of the asset and the end of its is the distance between Source i and MTF j and P i,j is the probability that a given casualty generated by Source i will be directed to MTF j (for clarity only the labels between Source 1 and MTF 1 are shown). participation in the simulation, and 2) a towed fuel level, which has a starting and maximum value of 400 fuel units and represents fuel available to be transferred to a site the asset visits.
Fuel levels in sites and assets can increase or decrease from one time step to the next. Fuel levels decrease through consumption of fuel, where each site consumes one fuel unit per time step and each asset consumes two fuel units per time step. Site fuel levels increase when a visiting asset transfers its towed fuel to the site and asset fuel levels increase when the asset visits a depot as described below.
The assets all begin the simulation in the depot. In each time step of the simulation, each asset performs one action. If the asset is in a depot and its internal fuel level is less than its maximum value, then the only possible action for the asset is to remain at the depot and increase its internal fuel level at a rate of two fuel units per time step. The asset's towed fuel value is set to its maximum value after the first time step the asset spends in the depot. If the asset is in a depot and its internal fuel level is at its maximum value, then the asset moves to a randomly determined neighboring site (site 1 or 31). If the asset is at a site, then the asset can perform one of up to three types of actions: 1) transfer some number of fuel units to the current site, where the number of fuel units transferred cannot exceed either the asset's current towed fuel value or be such that the site would exceed its maximum fuel value, 2) move from the current site to a neighboring site, or 3) stay at the current site and do not transfer any fuel. The asset randomly chooses one action from a list of possible actions formed by all legal varieties of the three action types, i.e., the asset has an equal probability of staying at the current site, transferring one unit of fuel, transferring two units of fuel, moving to one neighboring site, or moving to another neighboring site, assuming each action was legal for the particular circumstances. Because the scale of fuel values used in the simulations was in the tens or hundreds and the number of neighbors for each site was six or less, the expected asset behavior was to stay at a site transferring fuel for each time step until the site was near its fuel capacity.

Koopman Mode Analysis of MTF simulation results
For each value of N c (the casualty source rate), one hundred simulations were run, each with random MTF capacities subject to the constraints described in the previous section. Each simulation was run for 365 time steps (i.e. one year of simulated time). The observable value used as the KMD input was the MTF "fullness," defined as the ratio of the occupancy of each MTF to its capacity (see Fig 3 for example time series). The KMD algorithm used was the DMD_RRR method [23]. The output quantity of interest was the total Dead Of Wounds (DOW) for each simulation, where, clearly, a lower DOW was better. Fig 4 shows the mean and standard deviation for the 100 DOW values from each simulation for each casualty source rate case. The increase of the mean of DOW with casualty source rate is expected, as any MTF occupancy above the site's Table 3. Characteristic parameters of each MTF. P DOW is the probability that a given casualty in the MTF will die of wounds in a given time step, P RTD is the probability that a given casualty in the MTF will return to duty in the given time step, and the R out rate is the number of casualties (up to the total number of casualties present in the MTF) sent from the MTF to a hospital.   capacity is converted to DOW at the end of each day, so a larger casualty source rate naturally leads to a larger DOW. It is therefore interesting to instead consider a normalized DOW, defined as the mean DOW divided by the casualty source rate, which shows non-monotonic behavior. The standard deviation of the mean DOW shows the consistency of the resulting DOW for each simulation at a given casualty source rate, which is seen to approximately plateau at a high enough casualty source rate. In contrast, the standard deviation divided by the casualty rate shows non-monotonic behavior, as it has a peak around 130 casualties/source/ day, suggesting that the DOW value is most dependent on the relative MTF capacities at that casualty source rate and that MTF capacities need to be carefully determined to minimize the DOW value.
To gain greater insight into the change in system dynamics as the casualty source rate increases, it is instructive to apply Koopman mode decomposition of the fullness observable and examine the Koopman eigenvalue distributions. Fig 5 shows, for each casualty source rate, the distribution all of the Koopman eigenvalues from all 100 simulations at that casualty source rate. The eigenvalues show the expected behavior based on inspection of the MTF fullness time series. At low casualty source rate values, few or none of the MTF are near capacity so the time series show primarily damped oscillation around an equilibrium, which is reflected in the eigenvalue distribution by clustering around the origin. As the number of casualties from the sources increases, the MTF fullness values approach one and the damping effect decreases. Dynamically, this is a transition between a fluctuation-dissipation regime of fast decaying oscillations and a regime of slowly decaying or non-decaying fluctuations. As seen in the eigenvalue distributions, for increasing casualty source rate value, the imaginary component of the eigenvalues (corresponding to oscillatory behavior) tends toward the real axis, and the real component of the eigenvalues (corresponding to growth/decay behavior, and here showing the decrease in damping) tends toward more positive values.
The eigenvalue behavior is presented in a different form in Fig 6, which shows the mean and standard deviation of the real and imaginary components of the eigenvalues shown in Fig  5 and which show quantitatively the change in system dynamics described above. The mean real component is seen to increase to an approximately constant value with increasing casualty source rate, with a generally increasing standard deviation as well. The mean of the absolute value of the imaginary component and the standard deviation of the imaginary component are both seen to decrease with increasing casualty source rate.
Instead of considering all eigenvalues for each simulation, specific eigenvalues from each simulation and their corresponding modes can be selected. Two different eigenvalues/modes of interest are defined: the "dominant" mode and the "second" mode. These modes are defined as those whose corresponding eigenvalues have, respectively, the largest and second largest real components in each simulation.
The specific rational for the definitions of the dominant and second modes is that the longtime behavior of the site dynamics will by definition be dominated by those modes with the slowest decay rate (or equivalently, the fastest growth rate), as the time dependent coefficient of such modes in the Koopman decomposition will be larger than modes with smaller real components, over sufficiently long time scales. The mode defined as the dominant mode here will generally be the mode that is approximately equal to the mean of the data observables, and thus of limited dynamical interest. The mode defined as the second mode will better represent the mean-subtracted behavior of the system, as it includes a greater contribution of the damping effect on the system and oscillatory behavior. One could define further "higher order" modes using this approach, but the dynamical significance of such modes are less important on longer time scales, unless the real components of the associated eigenvalues are very similar to that of the second mode's eigenvalue, because of the exponential time dependence of the mode coefficients.
The dominant mode is generally expected to correspond to the mean occupancy of the MTFs, in which case it will be non-oscillatory and neither growing nor decaying. Figs 7 and 8 show the eigenvalues and their statistics of the dominant mode, which is indeed consistently at or near (1,0) and does not vary meaningfully based on the casualty source rate.
Analysis of the second mode is more useful to identify the qualitative change in behavior with increasing casualty rate. Figs 9 and 10 show the eigenvalues of the second mode and their statistics. A transition between the low casualty source rate and the high casualty source rate is apparent between casualty source rates of 30 to 110, corresponding to the transition in normalized DOW seen previously in Fig 4. This appears as the second mode eigenvalues transition from strongly damped oscillatory behavior to non-decaying fluctuations. Also the standard deviations of both the real and imaginary components decrease with increasing casualty source rate, corresponding to a transition to more similar behavior between random realizations with varied MTF capacities when the MTFs are consistently at or near capacity.
The power of KMD analysis in this application is thus shown by the relation between the transition in the normalized DOW data and the behavior of the Koopman eigenvalues, particularly the second mode eigenvalue. For the general MTF case with arbitrary parameter values, the Koopman eigenvalues in general and the second mode eigenvalues in particular can be used as a diagnostic "green light, yellow light, red light" test for the state of the system. If the eigenvalues are near the origin, then the system is likely in the unsaturated, damped oscillation, low DOW state and the system is behaving as desired ("green light"). If the eigenvalues are in the intermediate region between the origin and (1,0), then the system is likely in a transition state between the low and high DOW states and potential problems may appear ("yellow light"). And if the eigenvalues are at or near (1,0), then the system is likely in the undesired, non-decaying fluctuation, high DOW state ("red light"). Fig 11 shows a schematic demonstration of this concept for second mode eigenvalues. In a simulation case, this information could be used to end suboptimal simulations to saving computing time. In a real-world MTF or other facility case, this information could be useful for decision makers to reallocate resources or take other corrective actions to minimize deaths or other loses.  categorize the casualty rate 30 case as being in the green light regime, as the eigenvalues remain near the origin indicating a strong damping effect; casualty rates 50 and 70 as in the yellow light regime, as the eigenvalue distribution has moved away from the origin; and higher casualty rate cases as being in the red light regime, as most or all of the second mode eigenvalues are near to unique circle, indicated weakened damping. The mission life for each simulation was also recorded and processed. Fig 15 shows the mean and standard deviation of the mission life of each of the 100 simulations for each initial fuel value. Also shown is the difference between the mean mission life and the initial fuel value for each site, to better show the deviation of the mean mission life from a linear increase. The simulation ends when any site reaches zero fuel value and so the initial fuel value is the minimum possible mission life. This also explains the approximately linear increase in mean mission life times for each maximum site fuel value, as any site not refueled by an asset within the minimum possible mission life will cause the simulation to end. The standard deviation of    For the case where the eigenvalue distribution is clustered in the green region (near the origin), the system is in a highly damped, fluctuation-dissipation regime and is thus likely "safe" and unlikely to exhibit extreme behavior. For the case where a significant number of eigenvalues are in the yellow region, the system damping is decreasing and the system may be approach a transition to more extreme behavior. When the eigenvalue distribution is heavily in the red region, the system is expected to be displaying fluctuations with limited or no damping and therefore extreme and potentially dangerous or damaging behavior.

Koopman Mode Analysis of ship fueling logistics simulation results
https://doi.org/10.1371/journal.pone.0222023.g011 vector corresponds to a single site and so an element having large magnitude indicates that the corresponding site's dynamical behavior can be described, at least in part, by that mode's eigenvalue.
The analysis consisted of examination of the distributions of three sets of eigenvalues: 1) all of the eigenvalues from each simulation for each maximum site fuel value, 2) a single eigenvalue corresponding to the so-called "dominant" mode from each simulation for each maximum site fuel value (i.e. 100 eigenvalues per maximum site fuel value), and 3) a single eigenvalue corresponding to the so-called "second" mode. As for the MTF case, the dominant and second modes were defined, respectively, as the mode with the largest real eigenvalue and the mode with the second largest real eigenvalue, where these definitions were chosen to capture the dominant long-time behavior of the site dynamics.  Fig 23 shows that sites that are refueled one or more times by a significant amount (specifically, sites 1, 2, 3, 8, 9, 13, 15, 20, 21, 22, 26, 27, 28, 31, 32, 33 and 34) have relatively small positive or negative values in the dominant mode and relatively large values in the second mode, while the sites that have monotonically decreasing fuel values (sites 4, 5, 6, 7, 10, 11, 12, 14, 16, 17, 18, 19, 23, 24, 25, 29, 30, 35, and 36) have large positive values in the dominant mode and small positive values in the second mode. Because the dominant mode is a slowly decaying mode (i.e. its eigenvalue is just inside the unit circle) the larger positive mode values of the monotonically decreasing fuel sites correspond to primarily slowly decaying behavior, and the smaller positive or negative values of the refueled sites correspond to more slowly decaying or even growing behavior. Similarly, the second mode is further inside the unit circle and therefore represents faster decaying or more damped behavior, so the large values of the refueled sites represent damping of the applied forcing (i.e., the tendency of the fuel level to continue decreasing after being made to increase). Similar remarks apply to the 100 and 200 maximum site fuel value cases.
DMD decomposes signals into exponentially growing or decaying modes. However, in this problem, we have linearly decaying signals with random jumps. A priori, we have little intuition on what the true KMD spectrum is for such signals. As can be seen with the above results, DMD numerically computes a cluster of eigenvalues around 1 in the complex plane for such signals. Given this mismatch in the essential character of the signal (linear) and the approximating functions (exponential), one may wonder about the validity of the computed spectrum. In the next section, we perform a theoretical analysis of the algorithm that justifies the above numerical results. The analysis shows that using complex exponentials to approximate linear signals with random jumps forces repeated roots at 1. The clustering of eigenvalues around 1 is then due to numerical splitting of the repeated roots.

Approximating linear signals with KMD
In the ship fueling logistics problem, we have signals whose components decay linearly in time with random jumps in the signal. For these signals, we see a cluster of eigenvalues around 1 in the complex plane (see Fig 18). We argue here that this is due to a numerical splitting of repeated eigenvalues at z = 1. We show the analytic result for 2 random jumps. The argument can easily extended by induction. In what follows, we will be using notation from the Koopman Mode Decomposition section above. The simulation run times in this case are approximately twice as long as the 100 fuel unit case and consequently the assets have even more time to refuel sites, however it is seen that some site, such as 6, 12, 18, and 24, are still seldom or never refueled. Let 1 n ¼ ð1; . . . ; 1Þ 2 R n , and let c > 0 be the fuel capacity of the sites. Let t 1 < t 2 be positive integers which represent the random times when fuel is dropped at the sites. Let z 1 ; z 2 2 ðR þ Þ n be the amount of fuel dropped at the n sites at times t 1 and t 2 , respectively. The evolution of the fuel at the sites is given by where 0 < Δt � c is the time step of the simulation, H(t) is the Heaviside step function (=1 for t � 0 and 0 otherwise). We assume that m − 1 > t 2 .
In this case, X = (f 0 , . . ., f m−1 ) is spanned by The singular value decomposition of X can be expanded as where S is a 3 × 3 diagonal matrix, K 3 2 R n�3 , K nÀ 3 2 R n�nÀ 3 , W 3 2 R m�3 , and W mÀ 3 2 R m�mÀ 3 . The DMD matrix is given by where X 0 = (f 1 , . . ., f m ). Let us look more closely at the structure of A in light of (14). First, we recognize that the data is always spanned by B 3 = {1 n , z 1 , z 3 }, and the span of B 3 is the same as the subspace spanned by K 3 . Together, these imply that the image of U Δt is orthogonal to span K n−3 , and thus, 0 Expanding the definition of A: Therefore, the spectrum of A is Thus, we need to compute the spectrum of A 3 ¼ K � 3 U Dt K 3 to find σ(A). Analytic spectrum of A 3 . Note that A 3 is the representation of U Δt in the basis U 3 . Since the spectrum is basis-invariant, we can rewrite U Δt in terms of the basis B 3 . To do this, we use the relations • For the first basis element, we have Therefore, • For the second basis element, we have Furthermore, Equating (19) and (20) gives that • For the third basis element, we have Furthermore, ¼ ðc À t 2 DtÞ1 n À Dt1 n þ z 1 þ z 2 : ð24Þ Koopman Mode Analysis of agent-based models of logistics processes Equating the two expressions for U Dt f t 2 gives From (18), (21), and (25), we can write the matrix representation of U Δt in the B 3 basis as Since this is an upper triangular matrix, Finally, the total DMD spectrum is While this is the analytic spectrum, numerical errors due to floating point precision will give slightly different answers, as will be seen below.
Numerical example: DMD spectrum for linear decay. Here, we show that, although the analytically derived values for the DMD spectrum is given by (28)

Discussion and conclusions
The application of KMA to the two ABM systems described here is intended to demonstrate the utility of KMA in revealing information about the present and future behavior of the system that is not apparent simply from examination of the observables of the system. It is not intended as an optimization or modeling technique by itself, and the example ABM systems described in this paper were chosen for their dynamical features rather than their fidelity or relevance to real-world SCM or health care systems. Nonetheless, it is instructive to examine existing applications of ABM to optimization problems in SCM and health care systems, and to consider the possible contributions of KMA to these types of models.
There exists an extensive literature on applications of ABM to optimization problems (see Barbati et al for a recent literature review [26]), including applications of ABM-based simulation-optimization to problems in SCM and health care. KMA, while not an optimization technique itself, can be applied to optimization problems in various ways. For example, Oremland and Laubenbacher [27] describe a model reduction technique for ABM optimization, which could be aided by KMA by identifying the most significant temporal and spatial behaviors of the full model to produce and validate an appropriate reduced order model.
The SCM literature includes a large number of applications of ABM. The reader is referred to, e.g. Arvitrida's [28] review of ABM in SCM, specifically in a collaboration context, and to the extensive literature review of Fuller et al [29]. The common problem inherent to the application of ABM to SCM and other similarly complex systems is how to determine, in a general sense, the model parameters required to accurately reproduce existing systems and/or to produce an optimal system. This problem can be described as the determination, explicitly or implicitly, of a cost or reward function, where for a given state of the system the function outputs the expected optimal decision or action. The computational cost, for a given algorithm, of solving such an optimization problem depends on the state space over which the optimization  process occurs. KMA can potentially improve the effectiveness of solving such optimization problems, where representation of the system state by decomposition of the past history of the system observables into Koopman eigenvalues and modes describing the temporal and spatial dependencies of the system dynamics simplifies the optimization process. As a specific example, Fuller et al [29] describe an ABM using learning agents where an efficient model of the oil industry value chain is produced without a priori or expert knowledge of optimal policies. In  their approach, the agents learn an appropriate reward function early in the simulation then apply that learned reward function to choose actions in the rest of the simulation. The optimality of the model is therefore based on how well the learned reward function for a given state of the system matches the true reward function, which in turn is based on how effectively the learning algorithm uses the state space information of the system. The effectiveness of the learning algorithm can potentially be improved by use of the Koopman eigenvalues and mode as inputs to the learning algorithm, in addition to the current system state. Such an approach  includes the greater information about the system response by condensing the system dynamics into single values, rather than requiring the details of the past system behavior to be determined from the entire set of previous states of the system.
There also exists literature on applications of ABM or discrete event simulation (DSE) based optimization to health care resource planning, which could potentially benefit from KMA-derived insights. Cabrera et al [30], Weng et al [31], Wong et al [32], and Zeinali et al [33] describe ABM or DES based approaches to create decision support systems for hospital staffing and resource allocation to reduce patient wait times or optimize resource usage. Common features of such systems are the stochastic nature of the arrival times and numbers of incoming patients, and the varied duration of the patients' stay in the hospital. In our analysis of the MTF system, applying KMA to the time series of patient numbers can provide actionable information on when the available patient capacity of a facility is likely to be exceeded, thus giving warning of the need to re-allocate resources to meet patient care demands. Yousefi et al [34] describe an emergency department resource planning model wherein a metamodel is constructed based on an ABM simulation, and genetic algorithms and neural network based methods are used to train the metamodel. Such an approach could benefit from application of KMA to both the metamodel construction task, where a sensitivity analysis of the time series information could inform the essential features of the ABM simulation to include in the metamodel, and also the Koopman eigenvalues and modes can be used as features to train upon which to train the neural networks, where those KMA derived features each have a simple time dependence and thus generally simplify the task of training a neural network to recognize the essential system dynamics [35], [36].
For both the MTF and SF logistics systems, analysis of the Koopman eigenvalues is seen to provide useful information about the state of the system. In both cases, analysis of the Koopman eigenvalues shows the change in system dynamics from one regime to another as a control parameter is varied, specifically from a strongly damped, fluctuation-dissipation regime to a weakly or non-damped oscillatory regime. In the MTF system, this change in dynamics is due to the increasing number of incoming casualties overwhelming the MTFs and causing the system to shift to a dangerous, high death rate state. For the ship fueling logistics system, increasing site fuel capacity leads to a positive outcome as the refueling assets are more effective as sites can last longer between being refueled.
It is easy to see that, if a signal is a combination of an exponential signals, DMD will compute the correct spectrum. However, for linear signals with random jumps-such as in the ship fueling logistics problem-it is not a priori obvious that DMD will correctly compute the spectrum or, indeed, if the spectrum is meaningful. We took up the justification of our numerical results for the ship fueling logistics problem with a theoretical analysis of the DMD algorithm applied to linearly decaying signals with random jumps. Such signals force the DMD spectrum to have repeated eigenvalues at z = 1 in the complex plane, with the number of repeated eigenvalues equal to the number of random jumps in the signal. Numerical Koopman Mode Analysis of agent-based models of logistics processes inaccuracies due to floating point arithmetic split these repeated eigenvalues resulting the cluster of eigenvalues around 1 that were seen in our results.
Our results indicate that tracking of bifurcations in complex logistic systems is possible using Koopman Mode Analysis. We have studied two different logistics systems here, with two different types of dynamics and concluded that indicators of performance can be built based on selected Koopman Modes. It would be of great interest to extend these types of studies using only observational data from logistics systems.

Data management
All datasets generated by our models are housed in a Zenodo public repository and can be found at: http://doi.org/10.5281/zenodo.2567599. The datasets possess their own DOI number and can be cited as [37].