Koopman Mode Analysis of agent-based models of logistics processes


Authors: James Hogg aff001;  Maria Fonoberova aff001;  Igor Mezić aff001;  Ryan Mohr aff001
Authors place of work: Aimdyn, Inc., Santa Barbara, CA, United States of America aff001;  University of California Santa Barbara, Santa Barbara, CA, United States of America aff002
Published in the journal: PLoS ONE 14(9)
Category: Research Article
doi: 10.1371/journal.pone.0222023

Summary

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.

Keywords:

Physical sciences – Mathematics – Algebra – Linear algebra – Eigenvalues – Applied mathematics – Algorithms – Systems science – Agent-based modeling – Bifurcation theory – Materials science – Materials – Fuels – Chemistry – Chemical reactions – Decomposition – Engineering and technology – Energy and power – Fuels – Research and analysis methods – Simulation and modeling – Algorithms – Agent-based modeling – Computer and information sciences – Systems science – Agent-based modeling – Bifurcation theory – Biology and life sciences – Psychology – Behavior – Social sciences – Psychology – Behavior

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 [49] 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] studied a highly-detailed model of goods transport between a port and inland regions by different levels of agents on a dynamic graph. Binmad and Li [14] use an ABM approach to demonstrate the benefits of a forgiveness method to improve reputation systems in online commerce and social communities. Stavash et al [15] present an ABM in a war-gaming context to model the logistics of military operations.

Analysis of the results of ABM can benefit from data processing techniques that reveal details of the system behavior not readily apparent from visual examination of the model outputs, such as the nature of the time-dependent dynamics of a system and sensitivity to agent parameters. Koopman Mode Analysis (KMA) is a data processing technique that lends itself to such analyses. In this work, KMA was applied to identify transitions between different types of dynamics in numerical simulations of two different agent-based models: one simulating a large-scale, emergency medical treatment facility (MTF) system and one simulating a transport logistics situation which can considered as a simple supply chain management (SCM) problem.

Koopman Mode Analysis is a mathematical technique based on the spectral theory of the Koopman operator [1619]. The Koopman operator gives a linear representation of a complex, nonlinear or stochastic process without any loss of information. It does this by looking at the evolution of observables on the state space of the process rather than tracking the full state directly. This is advantageous since it is often impossible, intractable, or of no interest to fully specify the state space or the process—one merely needs to look at how observables on the system evolve in time.

Having its roots in harmonic analysis, KMA decomposes an observable’s signal into components having simple temporal dynamics plus a term dealing with the so-called continuous spectrum that one can think of as the noise term. It is analogous to how a Fourier transform decomposes a signal into sinusoids with specified frequencies and phases. KMA, however, gives more. First, in addition to identifying purely oscillating components like the Fourier transform, KMA can identify growing or decaying temporal signatures (with or without oscillation). Thus KMA components have a temporal signature given by a complex exponential rather than a sinusoid. Second, KMA has a spatial component. To each temporal signature, there is identified a complex spatial structure. Thus KMA represents an observable’s signal as a linear combination of these spatial structures (called “modes”) multiplied by their temporal coefficients. An important feature of KMA is that it can be deployed to data-only systems, namely those that do not have a model associated with them. However, in this paper we apply it to outputs of two ABM systems on which detailed understanding of the impact of the underlying methodology can be discerned.

The first ABM system analyzed in this paper consists of medical facilities receiving and outputting battlefield casualties, where the input rates of casualties and the patient capacity and treatment efficacy of the medical facilities determines the dynamics of patient survival rates. The second system involves mobile fuel delivery assets (representing, e.g., ships) traversing a network of fuel-consuming sites, the dynamics of which depend on the site fuel capacity and the ability of the assets to keep sites supplied with fuel.

The two systems displayed two different kinds of dynamics, one of which was primarily stationary with fluctuations about a mean, and the other dominated by linear decay. Linear decay dynamics have not previously been treated using KMA methods, and here we show a theoretical understanding of it in terms of degeneracy of eigenvalues. From these numerical and theoretical results, we show that KMA provides an effective method to study qualitative transitions in large logistics systems.

This paper is organized as follows: In the Materials and Methods section, we first provide an overview of Koopman Mode Analysis (KMA) and the Dynamic Mode Decomposition (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.

Materials and methods

Koopman Mode Decomposition (KMD)

Koopman Mode Analysis (KMA) decomposes a complex signal, { g t ∈ R n : t ∈ I ⊂ R }, 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 = aj + j is an eigenvalue of the underlying family of operators, a complex frequency determining the temporal signature; vj is the mode which determines the spatial structure of the associated temporal signature; the overlines denote complex conjugation; and n t ∈ 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 { g t ∈ R n : t ∈ I } is a stochastic process. There is an underlying state space Z, a flow St: ZZ (tI), and a vector-valued function g : Z → R n such that


The composition operation on the right-hand side induces a family of linear, infinite-dimensional operators (called the family of Koopman operators) U t : F → F, where F = { g : Z → R n } 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 {Ut} is a family of linear operators, we can analyze its spectrum. Assuming its a spectral family, the family can be decomposed as


Using this with (4), gives (1) (where we have suppressed the dependence of the vj’s and nt on g).

In KMA, the eigenvalues, e λ j t, and the modes, vj, from (4) are determined via the spectral analysis of a data-driven approximation of the family of Koopman operators Ut.

In terms of a logistics problem such as the Medical Treatment Facility (MTF) problem, the components of gt could represent the number of patients at each of the n medical facilities at time t. The state space and flow map (Z and St, respectively) would correspond to the state space for the underlying ABM model and its simulation, respectively. A Koopman mode v j = ( v j ( 1 ) , … , v j ( n ) ) T would determine which of medical facilities had its level of patients vary with a temporal signature partially characterized by the Koopman eigenvalue λj = aj + j. For example, if v j ( i ) was non-zero, then this would indicate that the number of patients at facility i is partially characterized by growth (aj > 0) or decay (aj < 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 t ( i ), v j ( i ), and n t ( i ) are the ith components of gt, vj and nt, 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, (f0, …, fm), where fk = gkΔt, as an input and returns the matrix V = (v0, …, vm−1) and the vector Γ = (γ0, …, γm−1). The i-th column of V and the i-th element of Γ are approximations of some vj from (1) and the associated eigenvalue e λ j Δ t, respectively. Note that fk = UΔt fk−1. The DMD algorithm first looks for a matrix A : R n → R n that approximates UΔt on the Krylov subspace K = span { f 0 , … , f m - 1 } ⊂ F. That is, we are looking for a matrix representation A of P K U Δ t : 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, vj and e λ j Δ t, respectively.

Algorithm 1: DMD

[V, Γ] = dmd(F)

Input: F = ( f 0 , … , f m ) ∈ R n × ( m + 1 )

1: X ← (f0, …, fm−1), X′ ← (f1, …, fm)

2: [K, Σ, W] = svd(X)

3: Determine numerical rank, k

4: KkK(:, 0 : k − 1), WkWk(:, 0 : k − 1), Σk ← Σ(0 : k − 1, 0 : k − 1)

5: A ← K k * X ′ W k Σ k - 1, (Schmid’s formula for K k * U Δ t K k in the SVD basis)

6: [Sk, Γ] = eig(A), (Γ = (γ0, …, γk−1), ASk(:, i) = γiSk(:, i), ‖Sk(:, i)‖2 = 1)

7: VKkSk

Return: V, Γ

Remark. Note, we are indexing from 0 in Algorithm 1, Kk is an orthogonal basis for the k-dimensional Krylov subspace K, and A is the matrix representation of P K U Δ t : K → K in the SVD basis Kk. The Koopman eigenfrequency λj can be computed from γj via λ j = 1 Δ t log ( γ j ), where we take the principal value of the complex logarithm.

KMA and logistics problems

In terms of the logistics problems, fk’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.

  • Measure quantities of interest, fk, at time tk (either from the simulation of an underlying ABM model or from a real system).

  • Input finite windows of snapshots (fkm, …, fk) into the DMD algorithm and return Koopman modes and eigenvalues vj and λj for that time window.

  • 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 changed—the 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 < hg < hr ≤ 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, hg, hr 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 Nc of casualties per day, where Nc = 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.

<h2>MTF network connectivity.</h2>
Fig. 1.

MTF network connectivity.


Di,j is the distance between Source i and MTF j and Pi,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).

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.

Tab. 1.

Distance (in km) between each source and each MTF.

<h2>Distance (in km) between each source and each MTF.</h2>
Tab. 2.

Probability that a casualty from a given source will be directed to a given MTF.

<h2>Probability that a casualty from a given source will be directed to a given MTF.</h2>
Tab. 3.

Characteristic parameters of each MTF.

<h2>Characteristic parameters of each MTF.</h2>

PDOW is the probability that a given casualty in the MTF will die of wounds in a given time step, PRTD is the probability that a given casualty in the MTF will return to duty in the given time step, and the Rout rate is the number of casualties (up to the total number of casualties present in the MTF) sent from the MTF to a hospital.

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.

<h2>Diagram of mission site locations and depot connectivity.</h2>
Fig. 2.

Diagram of mission site locations and depot connectivity.


The sites are on a 6-by-6 hexagonal grid and sites that share an edge are connected. In addition, mission sites 1 and 31 are connected to the depot as indicated by the arrows.

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 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.

Results

Koopman Mode Analysis of MTF simulation results

For each value of Nc (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].

<h2>Example MTF occupancy time series, for the 90 casualties / source / day case.</h2>
Fig. 3.

Example MTF occupancy time series, for the 90 casualties / source / day case.


The horizontal axis is days and the vertical axis is MTF occupancy. The vertical axis limit of each plot is one patient more than the capacity of that MTF, therefore the fullness at each time step is approximately the ratio of the plotted occupancy to the vertical limit of the plot. E.g., MTFs 2 and 10 have fullness values equal to or close to 1 for most of the simulated period, MTFs 9 and 13 often reach a fullness of 1, and the remaining MTFs have fullness values that varying around roughly 0.5.

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 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.

<h2>Mean and standard deviation of total Dead Of Wounds (DOW) for each casualty source rate.</h2>
Fig. 4.

Mean and standard deviation of total Dead Of Wounds (DOW) for each casualty source rate.


The top row shows the mean and standard deviation of DOW, and shows that the mean DOW increases superlinearly with casualty source rate, while the standard deviation reaches a maximum and levels out around a casualty source rate of around 175. To show the deviation from a linear increase, the bottom row shows the mean and standard deviation divided by the casualty source rate. This shows a sigmoidal shape to the mean plot and a clear peak to the standard deviation plot at around 130 casualties/source/day.

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.

<h2>All eigenvalues for all simulations for each casualty source rate.</h2>
Fig. 5.

All eigenvalues for all simulations for each casualty source rate.


The color scale (blue to red) indicates the L-2 norm of each eigenvalue’s associated Koopman mode. The transition from damped behavior at low casualty source rates to less damped behavior at higher casualty source rates is apparent in the shift in the eigenvalue distribution from clustering around the origin to collapsing onto the real axis and approaching the unit circle.

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.

<h2>Mean and standard deviation of the real and imaginary components of the eigenvalues in <em class="ref">Fig 5</em>.</h2>
Fig. 6.

Mean and standard deviation of the real and imaginary components of the eigenvalues in Fig 5.


The real components are seen to increase with casualty source rate and approach an approximate maximum, and the imaginary components are seen to approach the real axis.

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 long-time 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.

<h2>Eigenvalue distribution of dominant mode (i.e. mode with largest real eigenvalue) for each casualty source rate.</h2>
Fig. 7.

Eigenvalue distribution of dominant mode (i.e. mode with largest real eigenvalue) for each casualty source rate.


The color scale (blue to red) indicates the L-2 norm of each eigenvalue’s associated Koopman mode. The dominant mode eigenvalues are seen to not vary significantly with casualty source rate, consistent with their identification as the mean patient value over the MTFs.

<h2>Mean and standard deviation of eigenvalue distribution of dominant mode for each casualty source rate.</h2>
Fig. 8.

Mean and standard deviation of eigenvalue distribution of dominant mode for each casualty source rate.


The dominant mode eigenvalues are seen to not varying significantly with 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.

<h2>Eigenvalue distribution of second mode for each casualty source rate.</h2>
Fig. 9.

Eigenvalue distribution of second mode for each casualty source rate.


The color scale (blue to red) indicates the L-2 norm of each eigenvalue’s associated Koopman mode. The second mode eigenvalues are seen to represent the general dependence on casualty source rate of the all eigenvalue case shown in Fig 5.

<h2>Mean and standard deviation of eigenvalue distribution of second mode for each casualty source rate.</h2>
Fig. 10.

Mean and standard deviation of eigenvalue distribution of second mode for each casualty source rate.


As in the all-eigenvalue case, the real component approaches the unit circle as the imaginary component approaches the real axis.

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. Fig 11 shows a schematic illustration of this concept, where the boundaries between the colored regions are chosen arbitrarily. For a specific application the boundaries of the regions would be determined based on prior data and/or theoretical considerations. By way of example, in Fig 9, one might reasonably 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.

<h2>Schematic illustration of the possible use of the Koopman eigenvalue distribution as a diagnostic tool of the “health” of a dynamical system, where as noted in the text the specific radii of the colored regions shown here are arbitrary.</h2>
Fig. 11.

Schematic illustration of the possible use of the Koopman eigenvalue distribution as a diagnostic tool of the “health” of a dynamical system, where as noted in the text the specific radii of the colored regions shown here are arbitrary.


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

For each of the eight values of site maximum fuel capacity (25, 50, … 200), one hundred simulations were run and for each simulation the site fuel values at each time step were the observables used for Koopman Mode Analysis using a DMD-based algorithm (DMD_RRR). Figs 12, 13 and 14 show example site fuel value time series for single simulations with maximum site fuel values of 25, 100, and 200 fuel units. It is seen that the larger maximum site fuel values enable longer mission life times and more opportunities for assets to refuel sites, thus leading to more complicated dynamical behavior.

<h2>Time series of site fuel unit values, for 25 initial (and maximum) fuel units case.</h2>
Fig. 12.

Time series of site fuel unit values, for 25 initial (and maximum) fuel units case.


The linear decrease in fuel level is due to the consumption by each site of one fuel unit per time step in the simulation. The increases occasional increases in fuel level at some sites are due to refueling by mobile assets.

<h2>Time series of site fuel unit values, for 100 initial (and maximum) fuel units case.</h2>
Fig. 13.

Time series of site fuel unit values, for 100 initial (and maximum) fuel units case.


The longer simulation run times made possible by the higher initial site fuel value gives the assets more time to refuel sites.

<h2>Time series of site fuel unit values, for 200 initial (and maximum) fuel units case.</h2>
Fig. 14.

Time series of site fuel unit values, for 200 initial (and maximum) fuel units case.


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.

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 mission life times for each maximum site fuel value is also consistent with the increased possibility of the assets refueling all sites at least once, leading to mission life times longer than the minimum value. The deviation of the mean mission life values from the linear, minimum mission life values shows the same behavior as the mission life standard deviation, where larger maximum site fuel values increase the likelihood of all sites being refueled at least once before reaching zero fuel value.

<h2>Mission life statistics.</h2>
Fig. 15.

Mission life statistics.


The mean mission life increases with the maximum and initial site fuel value, and the standard deviation of mission lives also generally increases with maximum and initial site fuel value. To show the deviation from a purely linear increase in mission life, the bottom left plot shows the mean mission life with the minimum possible mission life subtracted off. This result shows the mission-life extending effect due to refueling by the assets.

Koopman Mode Analysis (KMA) processing was performed on the ship fueling logistics simulations to study the change in dynamical behavior as the maximum site fuel value was increased. A DMD-based algorithm for Koopman Mode Decomposition (KMD) was applied to the site fuel value time series from each of the 100 simulations for each maximum site fuel value, i.e. for each of the eight maximum site fuel values, KMD was performed on each of the 100 simulations, for a total of 800 applications of KMD. For each simulation, the KMD input consisted of a 2-D matrix with 36 columns and a number of rows each to the mission life in number of time steps. Each row consists of the fuel values of each of the 36 sites during a particular time step. The outputs of the KMD for each simulation are an equal number of Koopman modes and eigenvalues, where the mode is a 36 component vector and each eigenvalue is a single complex number that determines the time-dependence of the mode, specifically its real part determines the exponential growth/decay behavior and the imaginary component determines the oscillatory behavior of the mode. Examination of the magnitude of the elements of each mode gives insight into the behavior of each mode, as each element in a mode 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 16 shows all the Koopman eigenvalues for each maximum site fuel value and Fig 17 shows the mean and standard deviation of the real and imaginary components of each eigenvalue. These results show a transition from more strongly damped behavior at lower values of maximum site fuel value to less strongly damped behavior at higher values of maximum site fuel values.

<h2>All eigenvalues for each maximum site fuel value.</h2>
Fig. 16.

All eigenvalues for each maximum site fuel value.


The eigenvalue color scale (blue to red) shows the L-2 norm of each eigenvalue’s associated Koopman mode. Similar to the MTF case, the eigenvalues are seem to cluster toward the real axis and unit circle with increasing maximum site fuel values.

<h2>Mean and standard deviation of all eigenvalues for each maximum site fuel value.</h2>
Fig. 17.

Mean and standard deviation of all eigenvalues for each maximum site fuel value.


The real component is seen to go towards the unit circle and the imaginary component goes toward the real axis with increasing maximum site fuel values.

Fig 18 shows the eigenvalues for the dominant mode for each site initial fuel value and Fig 19 shows the mean and standard deviation values of the real and imaginary components of those eigenvalues. Figs 20 and 21 show the eigenvalues and their statistics for the second mode for each site initial fuel value. The eigenvalues of the dominant mode are seen to consistently cluster near (and generally inside) the unit circle with small or zero values of the eigenvalue imaginary component. These modes represent the constant (mean) or slowly growing/decaying components of the site fuel level dynamics. Additional information about the site dynamical behavior is apparent in the distribution of the second mode eigenvalues, where the transition from eigenvalues being generally closer to the origin to clustering near the dominant mode eigenvalues as the maximum site fuel value increases represents a decrease in damping or an increase in external driving. This is consistent with the increased occurrences of refueling made possible by the increasing minimum mission life, where the refueling function of the assets serves as a driving force.

<h2>Eigenvalues of the dominant mode for each simulation for each initial fuel value.</h2>
Fig. 18.

Eigenvalues of the dominant mode for each simulation for each initial fuel value.


The eigenvalue color scale (blue to red) shows the L-2 norm of each eigenvalue’s associated Koopman mode. The dominant mode eigenvalues mostly represent the mean of the site fuel values and do not depend strongly on maximum site fuel value.

<h2>Mean and standard deviation of dominant mode eigenvalue components.</h2>
Fig. 19.

Mean and standard deviation of dominant mode eigenvalue components.


The eigenvalues do not show a strong dependence on maximum site fuel value.

<h2>Eigenvalues of second mode for each simulation for each initial fuel value.</h2>
Fig. 20.

Eigenvalues of second mode for each simulation for each initial fuel value.


The eigenvalue color scale (blue to red) shows the L-2 norm of each eigenvalue’s associated Koopman mode. As in the MTF case, the second mode eigenvalues represent the general dependence of the all-eigenvalue case on the maximum site fuel value.

<h2>Mean and standard deviation of second mode eigenvalue components.</h2>
Fig. 21.

Mean and standard deviation of second mode eigenvalue components.


The second mode eigenvalues are seen to approach the unit circle and real axis with increasing values of maximum site fuel value.

Figs 22 to 27 show the dominant and second eigenvalues and corresponding modes (in bar graph and site spatial location form) for single example simulations for maximum site fuel values of 25, 100, and 200 units. Examination of the magnitudes of the mode elements reveals information about the dynamical behavior of each element’s corresponding site. For example, comparison of the time series in Fig 12 with the modes in 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.

<h2>Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 25 units.</h2>
Fig. 22.

Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 25 units.


The dominant mode eigenvalue is nearly on the unit circle and the real axis, indicating little or no temporal growth/decay or oscillatory behavior, while the second mode eigenvalue is closer to the origin and off of the real axis, indicating long term decay and oscillatory behavior of the associated mode.

<h2>Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 25 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).</h2>
Fig. 23.

Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 25 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).


The dominant mode values approximate the mean site fuel values, while the second mode values show the degree to which site fuel values oscillate and/or undergo long-term decay.

<h2>Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 100 units.</h2>
Fig. 24.

Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 100 units.


The dominant mode eigenvalue is nearly on the unit circle and the real axis, indicating little or no temporal growth/decay or oscillatory behavior, while the second mode eigenvalue is closer to the origin and off of the real axis, indicating long term decay and oscillatory behavior of the associated mode.

<h2>Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 100 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).</h2>
Fig. 25.

Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 100 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).


The dominant mode values approximate the mean site fuel values, while the second mode values show the degree to which site fuel values oscillate and/or undergo long-term decay.

<h2>Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 200 units.</h2>
Fig. 26.

Example eigenvalues of the dominant and second modes from a single simulation for maximum site fuel value of 200 units.


The dominant mode eigenvalue is nearly on the unit circle and the real axis, indicating little or no temporal growth/decay or oscillatory behavior, while the second mode eigenvalue is closer to the origin and off of the real axis, indicating long term decay and oscillatory behavior of the associated mode.

<h2>Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 200 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).</h2>
Fig. 27.

Example mode elements of the dominant and second modes from a single simulation for maximum site fuel value of 200 units, shown in bar graph form (a) and in the spatial positions of the corresponding sites (b).


The dominant mode values approximate the mean site fuel values, while the second mode values show the degree to which site fuel values oscillate and/or undergo long-term decay.

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.

Let 1 n = ( 1 , … , 1 ) ∈ R n, and let c > 0 be the fuel capacity of the sites. Let t1 < t2 be positive integers which represent the random times when fuel is dropped at the sites. Let z 1 , z 2 ∈ ( R + ) n be the amount of fuel dropped at the n sites at times t1 and t2, respectively. The evolution of the fuel at the sites is given by



where 0 < Δtc 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 > t2.

In this case, X = (f0, …, fm−1) is spanned by


The singular value decomposition of X can be expanded as


where S is a 3 × 3 diagonal matrix, K 3 ∈ R n × 3, K n - 3 ∈ R n × n - 3, W 3 ∈ R m × 3, and W m - 3 ∈ R m × m - 3. The DMD matrix is given by

where X′ = (f1, …, fm).

Let us look more closely at the structure of A in light of (14). First, we recognize that the data is always spanned by B3 = {1n, z1, z3}, and the span of B3 is the same as the subspace spanned by K3. Together, these imply that the image of UΔt is orthogonal to span Kn−3, and thus, 0 = K n - 3 * U Δ t X = K n - 3 * X ′. Expanding the definition of A:


Therefore, the spectrum of A is


Thus, we need to compute the spectrum of A 3 = K 3 * U Δ t K 3 to find σ(A).

Analytic spectrum of A3

Note that A3 is the representation of UΔt in the basis U3. Since the spectrum is basis-invariant, we can rewrite UΔt in terms of the basis B3. To do this, we use the relations UΔtf0 = f1, U Δ t f t 1 = f t 1 + 1, and U Δ t f t 2 = f t 2 + 1.

  • 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,


    Equating the two expressions for U Δ t f t 2 gives

From (18), (21), and (25), we can write the matrix representation of UΔt in the B3 basis as


Since this is an upper triangular matrix,


Finally, the total DMD spectrum is


For example, if Δt = 1 and the fuel capacity of the sites is c = 100, then the 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), the numerical computation gives slightly different values. For this example, we choose n = 50, m = 2n + 1, Δt = 1, and c = 100. The left panel of Fig 28 shows the first component of { f k } k = 0 m - 1, the middle panel the spectrum as computed by DMD, and the right panel shows the modulus of the spectrum. The numerical values of the spectrum of A3 computed using DMD_RRR are


Fig. 28.

(a) Linear decay with random jumps. (b) DMD_RRR spectrum corresponding to (30). (c) Modulus of the spectrum.

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 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].


Zdroje

1. Abar S, Theodoropoulos GK, Lemarinier P, O’Hare GMP. Agent Based Modelling and Simulation tools: A review of the state-of-art software. Computer Science Review. 2017;24:13–33. https://doi.org/10.1016/j.cosrev.2017.03.001.

2. Epstein JM, Axtell R. Growing artificial societies: social science from the bottom up. Brookings Institution Press; 1996.

3. Auyang SY. Foundations of complex-system theories: in economics, evolutionary biology, and statistical physics. Cambridge University Press; 1999.

4. Epstein JM. Modeling civil violence: An agent-based computational approach. Proceedings of the National Academy of Sciences. 2002;99(suppl 3):7243–7250.

5. Makowsky M. An agent-based model of mortality shocks, intergenerational effects, and urban crime. Journal of Artificial Societies and Social Simulation. 2006;9(2).

6. Rauhut H, Junker M. Punishment deters crime because humans are bounded in their strategic decision-making. Journal of Artificial Societies and Social Simulation. 2009;12(3):1.

7. Bosse T, Gerritsen C. Social Simulation and Analysis of the Dynamics of Criminal Hot Spots. Journal of Artificial Societies and Social Simulation. 2010;13.

8. Fonoberova M, Fonoberov VA, Mezic I, Mezic J, Brantingham PJ. Nonlinear dynamics of crime and violence in urban settings. Journal of Artificial Societies and Social Simulation. 2012;15(1):2.

9. Fonoberova M, Mezić I, Mezić J, Mohr R. An agent-based model of urban insurgence: Effect of gathering sites and Koopman mode analysis. PloS one. 2018;13(10):e0205259.

10. LeBaron B. Empirical regularities from interacting long-and short-memory investors in an agent-based stock market. Ieee transactions on evolutionary computation. 2001;5(5):442–455.

11. Huang CC. Using intelligent agents to manage fuzzy business processes. IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans. 2001;31(6):508–523.

12. Das R, Hanaoka S. An agent-based model for resource allocation during relief distribution. Journal of Humanitarian Logistics and Supply Chain Management. 2014.

13. Démare T, Bertelle C, Dutot A, Lévêque L. Modeling logistic systems with an agent-based model and dynamic graphs. Journal of Transport Geography. 2017;62:51–65. https://doi.org/10.1016/j.jtrangeo.2017.04.007.

14. Binmad R, Li M. Computational Models Based on Forgiveness Mechanism for Untrustworthy Agents. In: Lavangnananda K, Phon-Amnuaisuk S, Engchuan W, Chan JH, editors. Intelligent and Evolutionary Systems. Cham: Springer International Publishing; 2016. p. 29–42.

15. Stavash J, Chadha B, Wedgwood J, Welsh J, Parker M, Teitelbaum D. Agent Based Models for Logistics in Wargaming. In: Proceedings of the Fall 2003 SISO Simulation Interoperability Workshop; 2003.

16. Mezić I. Spectral Properties of Dynamical Systems, Model Reduction and Decompositions. Nonlinear Dynamics. 2005;41:309–325.

17. Mohr R, Mezić I. Construction of Eigenfunctions for Scalar-type Operators via Laplace Averages with Connections to the Koopman Operator. arXivorg. 2014; p. 1–25.

18. Mohr, R. M. Spectral Properties of the Koopman Operator in the Analysis of Nonstationary Dynamical Systems. Ph.D. thesis, University of California, Santa Barbara, Santa Barbara. 2014.

19. Korda M, Putinar M, Mezić I. Data-driven spectral analysis of the Koopman operator. arXiv. 2017;math.OC(1710.06532):1–18.

20. Schmid PJ, Sesterhenn J. Dynamic Mode Decomposition of Numerical and Experimental Data. In: Sixty-First Annual Meeting of the APS Division of Fluid Dynamics. San Antonio, Texas, USA; 2008.

21. Schmid PJ. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics. 2010;656:5–28.

22. Kutz JN, Brunton SL, Brunton BW, Proctor JL. Dynamic Mode Decomposition. Philadelphia, PA: Society for Industrial and Applied Mathematics; 2016.

23. Drmač Z, Mezić I, Mohr R. Data Driven Modal Decompositions: Analysis and Enhancements. SIAM Journal on Scientific Computing. 2018;40(4):A2253–A2285.

24. Mauroy A, Mezić I. Global stability analysis using the eigenfunctions of the Koopman operator. IEEE Transactions on Automatic Control. 2016;61(11):3356–3369.

25. Crnjaric-Zic N, Macesic S, Mezic I. Koopman Operator Spectrum for Random Dynamical Systems. arXiv preprint arXiv:171103146. 2017.

26. Barbati M, Bruno G, Genovese A. Applications of agent-based models for optimization problems: A literature review. Expert Systems with Applications. 2012;39(5):6020–6028. https://doi.org/10.1016/j.eswa.2011.12.015.

27. Oremland M, Laubenbacher R. Optimization of Agent-Based Models: Scaling Methods and Heuristic Algorithms. Journal of Artificial Societies and Social Simulation. 2014;17(2):6.

28. Arvitrida NI. A review of agent-based modeling approach in the supply chain collaboration context. IOP Conference Series: Materials Science and Engineering. 2018;337:012015.

29. Fuller D, Martins Ferreira Filho VJ, Arruda EF. Oil industry value chain simulation with learning agents. Computers & Chemical Engineering. 2018;111:199–209.

30. Cabrera E, Taboada M, Iglesias ML, Epelde F, Luque E. Optimization of Healthcare Emergency Departments by Agent-Based Simulation. In: ICCS; 2011.

31. Weng S, Tsai B, Wang L, Chang C, Gotcher D. Using simulation and Data Envelopment Analysis in optimal healthcare efficiency allocations. In: Proceedings of the 2011 Winter Simulation Conference (WSC); 2011. p. 1295–1305.

32. Wong SY, Tsui KL, Chin KS, Xu M. A simulation study to achieve healthcare service quality improvement in accident emergency department (AED). In: 2011 IEEE International Conference on Quality and Reliability; 2011. p. 259–263.

33. Zeinali F, Mahootchi M, Sepehri MM. Resource planning in the emergency departments: A simulation-based metamodeling approach. Simulation Modelling Practice and Theory. 2015;53:123–138. https://doi.org/10.1016/j.simpat.2015.02.002.

34. Yousefi M, Yousefi M, Ferreira RPM, Kim JH, Fogliatto FS. Chaotic genetic algorithm and Adaboost ensemble metamodeling approach for optimum resource planning in emergency departments. Artificial Intelligence in Medicine. 2018;84:23–33. https://doi.org/10.1016/j.artmed.2017.10.002.

35. Korda M, Mezić I. Learning Koopman eigenfunctions for prediction and control: the transient case. arXiv e-prints. 2018; p. arXiv:1810.08733.

36. Lusch B, Kutz JN, Brunton SL. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications. 2018;9.

37. Hogg J, Fonoberova M, Mezić I., Mohr R. Simulation results of two agent-based models of logistics systems. 2019. Available from: http://doi.org/10.5281/zenodo.2567599.


Článek vyšel v časopise

PLOS One


2019 Číslo 9