extfieldsofstudy:
Computer Science
Mathematics

text:
Efficient Capacity Computation and Power Optimization for Relay Networks

The capacity or approximations to capacity of various single-source single-destination relay network models has been characterized in terms of the cut-set upper bound. In principle, a direct computation of this bound requires evaluating the cut capacity over exponentially many cuts. We show that the minimum cut capacity of a relay network under some special assumptions can be cast as a minimization of a submodular function, and as a result, can be computed efficiently. We use this result to show that the capacity, or an approximation to the capacity within a constant gap for the Gaussian, wireless erasure, and Avestimehr-Diggavi-Tse deterministic relay network models can be computed in polynomial time. We present some empirical results showing that computing constant-gap approximations to the capacity of Gaussian relay networks with around 300 nodes can be done in order of minutes. For Gaussian networks, cut-set capacities are also functions of the powers assigned to the nodes. We consider a family of power optimization problems and show that they can be solved in polynomial time. In particular, we show that the minimization of the sum of powers assigned to the nodes subject to a minimum rate constraint (measured in terms of cut-set bounds) can be computed in polynomial time. We propose an heuristic algorithm to solve this problem and measure its performance through simulations on random Gaussian networks. We observe that in the optimal allocations most of the power is assigned to a small subset of relays, which suggests that network simplification may be possible without excessive performance degradation.

is the Gaussian relay network, which models power limited transmitters and received signals corrupted by additive white Gaussian noise.
While the capacity of some network models (e.g. wireless erasure and ADT) is well characterized, the capacity of the Gaussian relay network, even in its simplest form with one transmitter, one relay, and one receiver, is in general unknown. The best known capacity upper bound is the so-called cut-set bound. A cut Ω of a network can be considered as a subset of nodes which includes the source node and excludes the destination node. For this cut, the capacity F (Ω) is defined as the maximum rate that information can be transferred form the nodes in Ω to the nodes that are not in Ω conditioned on the fact the information on Ω c (the nodes that are not in Ω) is known. The cut-set upper bound is the minimum cut capacity over all the possible cuts.
In the Gaussian setting, there are several capacity lower bounds based on different communication schemes, such as amplify-and-forward, decode-and-forward, compress-andforward, quantize-and-forward, etc. [5], [7], [18]. Recently, Avestimehr, et al. [2] made significant progress in the capacity characterization of Gaussian relay networks by showing that a quantization and coding communication scheme can achieve a communication rate within a constant gap of the cut-set upper bound, where the gap only depends on the number of nodes in the network (i.e. it is independent of the channel gains and power levels). However, the evaluation of the achievable communication rate, which is necessary to implement the scheme, requires the computation of the cut-set bound for the network. Assuming that for a given cut the cut capacity is easy to compute, finding the cut-set upper bound can be a challenging problem. For a network with n relays there are 2 n different cuts and a greedy algorithm needs exponential time in the number of relays to compute the cut-set bound.
In this work we show that the achievable rate of the scheme of [2] for the Gaussian relay network can be computed in polynomial time, and as a result, can be computed efficiently. This result is obtained by showing that the cut capacity of a fairly large class of networks under the assumption of independent encoding at the nodes in Ω is a submodular function. For the special case of layered relay networks, [27] showed the equivalent of our submodularity result simultaneously with our conference version of this paper [25]. Submodularity properties of conditional entropy (in terms of which cutcapacities are expressed) have also been used in [2], [3] to bound the cut-capacity of a network in terms of the cutcapacity of the corresponding unfolded graph 1 .
Existing results on minimization of submodular functions provide algorithms with polynomial running time O(n 5 α + n 6 ), where α is the time that it takes to compute F (Ω) and n is the number of nodes in the network [20]. In addition, there exist possibly faster algorithms without polynomial time performance guarantees based on Wolfe's minimization norm algorithm [12]. In Section VI, by simulations, we show that the cut-set bound for a Gaussian relay network with around 300 nodes can be computed on a laptop computer in about a minute using a Matlab package for submodular minimization provided in [19].
Our results, extend and generalize previous results for the ADT model. This model can be seen as a high signalto-noise-ratio (SNR) approximation of the Gaussian model, incorporating the effects of broadcasting and superposition of signals while de-emphasizing the effects of noise. Amaudruz et al. [1] showed that the cut-set bound for a layered 2 ADT model can be computed efficiently. They have extended graph flow algorithms such as Ford-Fulkerson's in a nontrivial way to find the maximum possible linearly independent (LI) paths in the network. They showed that the capacity of the network is equal to the maximum number of (LI) paths and can be computed in time O(M · |E| · C 5 ), where M is the maximum number of nodes per layer, |E| is the total number of edges and C is the capacity of the network. Moreover, they showed that the capacity can be achieved by using a very simple one-bit processing at the relay nodes. Later Goemans et al. [13] showed that the deterministic model is a special case of a flow model based on linking systems, a combinatorial structure with a tight connection to matroids. As a by-product, they obtained the submodularity of the cut capacity for layered ADT networks. Using this observation they provided various algorithms related to matroid theory to compute the cut capacity of the layered deterministic model based on finding intersection or partition of matroids. These results led to faster algorithms to compute the capacity of large layered ADT networks. In addition, there has been other extensions on improving the running time of the current algorithms for computing the capacity of ADT networks [9], [10], [22], [23].
In addition to showing that the capacity within a constant gap of the Gaussian relay network can be computed in polynomial time, our results allow us to compute in polynomial time the capacity of the wireless erasure network. Furthermore, we provide a simple proof for the computability in polynomial time of the capacity of the layered and non-layered ADT networks.
Building on the submodularity of the cut-capacity for independent encoding at the nodes, we show that, in the Gaussian setting, it is possible to efficiently optimize the power allocated to the source and relay nodes. We consider two power optimization problems: (i) minimize the total power satisfying a minimum source-destination data rate constraint and power constraints at each node; (ii) maximize the sourcedestination data rate satisfying total and individual power constraints at the nodes. Since the capacity of the Gaussian relay network is approximately given by the cut-set upper bound with independent encoding at the nodes, we use this cut-set bound to characterize data rate in the optimization problems. We show that these optimization problems can be solved in polynomial time and use simulations to get insights about some properties of the optimal power allocations for networks of various sizes. We observe that optimal power allocations assign most of the power to a small subset of nodes and that setting the power to zero in the remaining nodes (i.e. removing these nodes from the network) often results in a small rate loss. Nazaroglu, et al. showed in [26] that for the special case of the N -relay Gaussian diamond network a fraction k/(k + 1) of the total capacity can be approximately achieved by using only k of the total N relays. This suggests that the diamond network can be significantly simplified by tolerating a small performance loss. Our results provide a numerical counterpart to the fundamental performance bounds derived in [26] and suggest that network simplification may also be possible in more general Gaussian relay networks.
We obtain these results by considering a general framework to compute the cut-set bound. We assign transmit signal random variable X i to node i ∈ V and we assume the probability distribution over the signals X 1 , X 2 , . . . , X n to be independent, i.e p(X 1 , X 2 , . . . , X n ) = p 1 (X 1 )p 2 (X 2 ) · · · p n (X n ). We also assign received signal random variables Y i 's to each node. The network is defined by the transition probability function f (Y 1 , Y 2 , . . . , Y n |X 1 , X 2 , . . . , X n ). We further assume that the transition probability function is of the form f 1 (Y 1 |X 1 , . . . , X n ) · · · f n (Y n |X 1 , . . . , X n ), meaning that the received signals are independent conditioned on the transmitted signals in the network. For such networks we show that F (Ω) = I(Y Ω c ; X Ω |X Ω c ) 3 is submodular with respect to Ω. Later we show that for ADT networks, the Gaussian relay network and the wireless erasure network, we can find p 1 (X 1 ) · · · p n (X n ) such that min Ω F (Ω) becomes equal to the capacity or the capacity within a constant gap. In other words, the min-cut problem for these networks can be cast as a minimization of a submodular function.
The paper is organized as follows. In Section III we show that for specific type of networks the cut value, F (Ω), is a submodular function. We then show in Section IV that for many wireless network models such as the ADT deterministic network, Gaussian relay network and wireless erasure network the capacity or an approximation to the capacity can be cast as a minimization of F (Ω). In Section V we study two power optimization problems and show that they can be solved efficiently. Finally, in Section VI we describe results related to solving optimization problems involving submodular functions and perform power optimization in various randomly generated networks of different sizes. We start by introducing the notation used in the rest of the paper.

II. NOTATION
Let V denote the set of nodes in the network and |V| its cardinality. For any subset A of nodes we denote by V\A or A c the set of nodes in V that are not in A. We assume V\A ∪ B = V\(A ∪ B). A cut Ω is defined as a subset of nodes in V. A cut splits the nodes in the network into two groups, the nodes that are in Ω and the ones that belong to V\Ω. Random variables are shown in capital letters such as X i and Y i . We use boldface letter for vectors, e.g. x is a constant vector and X is a random vector. We use X Ω to denote (X v1 , X v2 , . . . , X v |Ω| ) with v i ∈ Ω. The function I(X; Y |Z) is the mutual information between random variables X and Y conditioned on random variable Z. With a slight abuse of notation we use H(X) to denote either the entropy or differential entropy of the discrete or continuous random variable X [6]. By F p we denote a finite field with p elements. Finally, all the log(·) functions are in base two.

III. SUBMODULARITY OF CUT-SET FUNCTION
Submodularity arises in many combinatorial optimization problems and large body of research has been developed on minimizing or maximizing submodular functions under various constraints.
A submodular function f : 2 V → R is defined as a function over subsets of V with diminishing marginal returns, i.e. if A, B ⊆ V with A ⊆ B and any v ∈ V\B, The theorem below establishes the submodularity of the cut capacity function of a general relay network under some special assumptions. This theorem will be used in Section IV to prove that the capacity or an approximation to the capacity of various specific relay network models can be computed by minimizing a submodular function.

Theorem 1. Consider a network consisting of nodes in V.
Each node sends a message X i , i ∈ V and receives Y i , i ∈ V. If the messages are independent p(X 1 , X 2 , . . . , X |V| ) = p 1 (X 1 )p 2 (X 2 ) · · · p |V| (X |V| ) and conditioned on the sent messages the received messages are independent, then the function where (a) is the definition of mutual information and (b) is from the chain rule for the entropy function. Therefore, In the following example we show that if the signals at the nodes are correlated then F (A) is not necessarily a submodular function. Example. Consider a symmetric Gaussian diamond network with two relays such that the channel gains from source to relays are equal to one and from relays to destination are equal to three. Letting X s , X r1 , and X r2 be the signals transmitted at the source and relay nodes, then the received signals at relays and destination are given by 1). For this example, we set the probability distribution of X s , X r1 , X r2 to be jointly Gaussian with zero mean and covariance matrix Finally, consider the sets A = {s, r 1 } and B = {s, r 2 }. Figure 1 shows how the function depending on the value of ρ. It follows that in general F (·) is not a submodular or a supermodular function when the there is correlation among the signals at the nodes.

IV. WIRELESS NETWORK MODELS
In this section, by applying the result of Theorem 1, we show that the capacity or an approximation to the capacity for the ADT deterministic network, Gaussian relay network, and wireless erasure network can be cast as a minimization of a submodular function.

A. Deterministic model (ADT)
We start by briefly describing the network model of [4] 4 . In this model, each link from node i to node j has an associated non-negative integer gain n ij . Each node i ∈ V transmits a signal X i and receives a signal Y i , both in F q p where q = max i,j n ij . At any given time, the received signal at node j is given by where d is the destination node, the shifting matrix S is given by and the sums and products are in F p . For a given cut Ω of the network, where Ω includes the source node and excludes the destination node, we can stack together the input vectors X i , i ∈ Ω and output vectors Y i , i ∈ Ω c , and define a transition matrix Λ Ω that gives the input-output relationship of these vectors according to (1). It is shown in [2] that the capacity of the deterministic network is equal to min Ω rank(Λ Ω ). We show next in Theorem 2 that rank(Λ Ω ) is submodular, and hence the capacity can be computed by minimizing a submodular function.
, and let G be the set of cosets of N in F n p . Pickx i to be an element in the ith coset of N for i = 1, 2, . . . , |G |, and set y i = Ax i . Notice that y i = y j if i = j. Now, if we choose x uniformly at random from elements of F n p with probability 1/|F n p |, then the mapping Ax maps x to {y 1 , y 2 , . . . , y |G | } uniformly at random with probability 1/|G |. In addition, the cosets of N form a partition of F n p into p n /|N | sets. Also rank(A) + log p (|N |) = n. Thus, log p |G | = rank(A).

Theorem 2. For a deterministic model, given a cut Ω assume
Λ Ω is the transition matrix form nodes in Ω to nodes in Ω c . Set D(Ω) = rank(Λ Ω ), then D(Ω) is submodular.

Remark 1.
A special case of Theorem 2 for layered ADT networks was proved in earlier works [13], [23].
From the definition of transition matrix, Λ Ω , if we assume for the cut Ω, s = (s 1 , s 2 , . . . , s k ) t symbols are being sent from nodes in Ω and r = (r 1 , r 2 , . . . , r ℓ ) t symbols are being received by nodes in Ω c then r = Λ Ω s. Then we can write is the result of Proposition 1 and the fact the s has uniform probability distribution.
Notice that for the independent probability distribution on the sources the received signals are independent conditioned on transmitted signals so, based on Theorem 1,

B. Gaussian relay network
The Gaussian network model captures the effects of broadcasting, superposition and noise of power constrained wireless networks. In this model, at any time index (which we omit) the received signal at node j ∈ V\{s} is given by where X i ∈ C is the transmitted signal at node i, subject to an average power constraint E(|X i | 2 ) 1, h ij ∈ C is the channel gain from node i to node j, and N j ∈ CN (0, 1) is additive white circularly symmetric complex Gaussian noise, independent for different j.
It has been show in [21, Theorem 2.1] that using lattice codes for transmission and quantization at the relays, all rates R between source {s} and destination {d} satisfying can be achieved, where Ω is a source-destination cut of the network and X Ω = {X i , i ∈ Ω} are i.i.d . CN (0, 1). In addition, the restriction to i.i.d. Gaussian input distributions is within |V| bits/s/Hz of the cut-set upper bound [2]. Therefore the rate achieved using lattice codes in the above result is within 2|V| bits/s/Hz of the capacity of the network. The following corollary is an immediate consequence of Theorem 1.
Due to Corollary 1 the minimization in (3) is the minimization of a submodular function and the resulting optimal value is within 2|V| of the capacity of the network. 5 .

C. Wireless erasure network
In [8] the authors introduce a special class of wireless networks, called wireless erasure networks. In these networks, a directed graph G = (V, E) defines the interconnections between nodes. To model the broadcast effect of wireless networks, the signals on all outgoing arcs of any given node are equal to each other. There is no interference among multiple arcs arriving at a given node in this model, and the signals on the various arcs are erased independently of each other. We assume binary transmitted signals at each node, i.e. X i ∈ {0, 1}, i ∈ V\{d}, but all the results can be extended to models with larger input alphabets. It has been shown in [8] that the capacity of the network is where ǫ ij is the probability of erasure when node i is sending information to node j. We show in the following theorem that F (Ω) is submodular.
5 Notice that I(X Ω ; Y Ω c |X Ω c ) = log det(I + HH † ) where H is the matrix of channel gains from nodes in Ω to nodes in Ω c and H † is the conjugate transpose of H. Therefore, it is easy to compute the capacity of each cut.
We used in (a) the independence among X i and the channel erasures, in (b) the fact that for X i ∼ Bernoulli(1/2), H(X i ) = 1, and in (c) the fact that for X i ∼ Bernoulli(1/2), H(X i |Y j = e, j ∈ Ω c ) = 1 and for independent erasures we have p(Y j = e, j ∈ Ω c ) = j∈Ω c ǫ ij . Theorem 1 can be applied to conclude that F (Ω) is submodular.

V. POWER OPTIMIZATION
In the previous section, for the Gaussian relay network model, we considered fixed power assignments to the different nodes in the network, and have shown that a constant gap approximation to the capacity can be efficiently computed. In many applications it is of interest to allocate the nodes' transmission powers to optimize a given objective. For example, in a network where the nodes are battery powered, it may be of interest to maximize the network lifetime while satisfying a baseline quality of service. Alternatively, it may be desirable to maximize the network throughput for a given total power budget. This total power budget may arise due to, e.g., a maximum system weight constraint which is dominated by the battery weight, or a total system cost, which may be heavily influenced by the cost of the batteries. Power allocation optimization may also naturally arise in situations where the channel gains, while known to all the nodes, slowly vary over time. In this case, it may be desirable to optimally allocate power for the current channel condition.
As before, we characterize communication rates in terms of cut-set capacities. We consider a model where the cut-set capacities are functions of the cuts and powers assigned to the nodes in the network: F (Ω, p) : V × R |V| → R, and we focus on the Gaussian model where this function depends explicitly on the power assignment, F (Ω, p) = I(X Ω ; Y Ω c |X Ω c ) = log det(I + H Ω P Ω H † Ω ) (5) where H Ω is the matrix of channel gains from nodes in Ω to nodes in Ω c , H † is the conjugate transpose of H, and P Ω is a diagonal matrix where the diagonal elements are the powers of the nodes in Ω.
We will show in Lemma 1 below that in the Gaussian case F (Ω, p) is a concave function of p. While the results of this section are stated and proved for the Gaussian model, we conjecture that similar results should hold for other models in which F (Ω, p) is a concave function of p.
For Gaussian relay networks, we show that the following optimization problem can be solved in polynomial time, minimize R,P,p for fixed constants µ 1 , µ 2 , R 0 , p max and P tot . In the rest of the section, we denote the feasible set of the optimization (6) by K.
We use the Ellipsoid method [14], [15] to show that the optimization (6) can be solved efficiently. We will use the following definitions and result. The reader is referred to [15] for more details.

Definition 1 (Polynomial computability). A family of optimization programs is polynomially computable if:
(i) for any instance of the program and any point x in the domain, the objective and its subgradient can be computed in polynomial time in the size of the instance. (ii) for a given measure of infeasibility Infeas(·), it should be possible to determine if Infeas(x) ε in polynomial time, and when this inequality is not satisfied, it should be possible to find in polynomial time a vector c such that c T x > c T y, ∀y : Infeas(y) ε.

Definition 2 (Polynomial growth). A family of optimization programs has polynomial growth if the objectives and the infeasibility measures as functions of points
x in the domain grow polynomially with x 1 .

Proposition 2 ([15, Theorem 5.3.1]). Let P be a family of convex optimization programs equipped with infeasibility measure Infeas(·). Assume that the family is polynomially computable with polynomial growth and with polynomially bounded feasible sets. Then P is polynomially solvable.
In order to use Proposition 2 we need to check that the optimization (6) is a convex program. Since the objective function is linear, we only need to check that the feasible set K is convex.
It is easy to check that the sets P 1 = {(R, P, p) : 0 p p max }, P 2 = {(R, P, p) : p i P }, P 3 = {(R, P, p) : R 0 R}, and P 4 = {(R, P, p) : P P tot } are also convex sets, and, as a result, Having proved that (6) is a convex program, in order to use Proposition 2 we need to check that the conditions of Definitions 1, 2, and 3 are satisfied. Part (i) of Definition 1 follows from the linearity of the objective in (6). For part (ii) of Definition 1 we specify an infeasibility measure Infeas(·) : R |V|+2 → R as follows 6 : The conditions of part (ii) of Definition 1 are verified in the following theorem.
Proof: Part (a) requires checking that each of the arguments of the max of (7) is smaller than or equal to ε in polynomial time. The first six terms are linear functions and can be easily computed. The last term can be compared to ε by performing a minimization of a submodular function, which as was shown in Section IV-B, can also be computed in polynomial time.
We focus on condition (b). In this case Infeas((R, P, p)) > ε, meaning that at least one of the arguments of the max of (7) is larger than ε. We consider each case separately.
If −p i > ε we set c = −e i+2 where e i has a one in the i th position and zeros everywhere else, which can be easily checked to satisfy the condition of part (b).
Having proved these preliminary results, we are ready to prove the main result of this section.
Theorem 5. The optimization in (6) can be solved in polynomial time on the size of the problem.
Proof: The proof uses Proposition 2, which requires verifying the convexity of the problem together with the conditions of polynomial computability, polynomial growth, and polynomial boundedness of the feasible set. Convexity was proved in Lemma 1, while polynomial computability was shown in Theorem 4. Polynomial growth follows from the fact that F (Ω, p) is the log of a polynomial on p with degree at most |V|, while the objective and remaining terms that define the infeasibility measure are linear on (R, P, p). Finally, to check that feasible set if polynomially bounded, we note that the feasible set is a subset of the hypercube . It follows that the feasible set is contained in the Euclidean ball centered at the origin with radius (R max , P tot , p max ) 2 , which can be easily checked to grow polynomially on the size of the problem.
The general optimization problem (6) can be specialized to a power minimization with a minimum rate constraint, and to a rate maximization with a total power constraint. Both problems can be solved in polynomial time, as stated in the following corollaries to Theorem 5.

Corollary 2.
The following power minimization problem can be solved in polynomial time.

VI. ALGORITHMS AND SIMULATIONS
In this section we study different algorithms and and provide simulation results regarding submodular function minimization and power allocation problems. In the first subsection we look at minimum norm algorithm for submodular function minimization and we show that this algorithm can find the approximation to capacity of layered Gaussian relay networks with more that 300 nodes in a couple of minutes. In the second subsection we propose a heuristic algorithm to find the optimum power allocation for Gaussian relay networks.

A. submodular function minimization
One approach to solve the submodular minimization problem due to Lovász is based on extension of the set function f : 2 V → R to a convex function g : [0, 1] |V| → R that agrees with f on the vertices of the hypercube [0, 1] |V| , with a guarantee that min A⊆V f (A) is equal to min x g(x) for x ∈ [0, 1] |V| . In this section we assume the normalization f (∅) = 0.
The Lovász extension g of any set function f can be defined as follows. For a given x ∈ [0, 1] |V| order the elements of A key result is that if f is submodular its Lovász extension g is a convex function [14], [11]. In addition, finding the minimum of the submodular function f over subsets of V is equivalent to finding the minimum of the convex function g in the hypercube [0, 1] |V| . The optimization can be done in polynomial time using Ellipsoid algorithm [14].
There are other algorithms with faster running time to solve the submodular minimization problem [16], [17], [20]. To the best of our knowledge, the running time of the fastest algorithm is in the order of O(n 5 α + n 6 ), where α is the time that the algorithms takes to compute f (A) for any subset A ⊆ V [20]. For ADT networks, Gaussian relay networks, and erasure networks, α is the time to compute: the rank of n × n matrices, the determinant of n × n matrices, and equation (4), respectively.
However, for networks of large size, a complexity of O(n 5 α + n 6 ) may still be computationally cumbersome. As a result, in these cases it is desirable to have faster algorithms. Recently, Fujishing [11], [12] showed that the minimization of any submodular function can be cast as a and the corresponding minimum norm optimization is Letting x * be the solution of this minimization, the set A * = {v i : x * (v i ) < 0} is the solution to min A f (A). Whether the above optimization problem can be solved in polynomial time is an open problem. However empirical studies [12] have shown that this algorithm has comparable or even faster running times than the other algorithms with polynomial time performance guarantees. In our specific setting, for layered Gaussian relay networks of size up to around 300 nodes with 4 nodes per layer, we were able to find the approximate capacity (cf. (3)) in order of minutes (see Figure 2). In order to solve the minimization (10) we used the Matlab package provided in [19].

B. Power allocation
In Section V we have shown that the Ellipsoid method can be used to solve the optimization in (8) in polynomial time.
While in theory this result shows that the optimization in (8) is tractable, in practice the Ellipsoid method has a number of shortcomings that limit its usability. On the one hand the running time of the Ellipsoid method can be large compared to alternative algorithms, and on the other hand, for high dimensional problems it has shown numerical instability. In this section, we propose a heuristic algorithm to solve the optimization in (8) and show that this algorithm converges to the right solution. While the running time can be exponential on the network size, we show through simulations that the algorithm often converges within a time proportional to the network size.
Our proposed algorithm, Algorithm 1 (see pseudo-code below), is based on the cutting plane methods [24] in convex optimization. The optimization in (8) contains exponentially many constraints of the form In Algorithm 1 we first find the min-cut corresponding to assigning maximum power to all nodes in the network: Then, we modify the optimization in (8) by replacing the constraint (11) with The resulting convex program can be easily solved since it contains few constraints. After optimization, we let p * be the optimum power allocation for the current set of constraints and set We iteratively add the constraint to our set of constraints, and solve the optimization again. We stop if the new constraint is already in the set of constraints.

Algorithm 1 Power minimization
Input: Channel gain matrix H, desired rate R, vector of nodes' power constraints p max . Output: Min-cut, Power assignment p * that achieves approximate to rate R with minimum sum of powers.
Since in each iteration the algorithm adds a new constraint to the constraint set and the number of constraints in (11) is finite, the algorithm is guaranteed to find the optimum power in a finite number of iterations, which can be exponential.
We used simulations to test the performance of the algorithm for networks of varying size n ranging from 10 to 40. For each n, we generated 300 random networks with channels gains drawn i.i.d. using a N (0, 1) distribution. We set the desired transmission rate R 0 = 4 in (8) and set a maximum power constraint in each node to p max = 100. The results are  shown in Figures 3-6, where the the vertical bars represent ±1 standard deviation around the mean computed over the 300 random networks. Figure 3 shows the number of iterations of Algorithm 1 as a function of the number of nodes in the network. We see that the number of iterations grows O(n 3/2 ) with the number of nodes.
In Figure 4, we present the simulation time as a function of the network size. We observe a running time that grows slower than O(n 6 ). The figure also shows that the power optimization of networks of 40 nodes completes in less than two hours on an Intel Xeon quad core CPU running at 2.33 GHz.
The minimum sum of powers to approximately achieve R 0 = 4 for networks of different size is presented in Figure 5 with a blue solid line. Interestingly, the plot shows that the minimum power concentrates around the mean. In addition, the figure shows diminishing returns in total power savings resulting from increasing the network size.   For the special case of the diamond relay network with N relays, [26] shows that a fraction k/(k + 1) of the network capacity can be approximately achieved by using only k relays. In our setting of a general Gaussian relay network and sumpower minimization, we are interested in investigating whether it is possible to remove a large fraction of the nodes from the network without significantly affecting its performance.
In order to determine which nodes to remove from the network, we solve the sum-power minimization and compare the optimal power allocation p * i of each node i to a threshold P th . All relays with p * i < P th are removed from the network. Let N be the set of nodes with p * i P th . We optimize the power allocation for the network with node set {s} ∪ N ∪ {d} and determine whether the problem is feasible for R 0 = 4. If the problem is infeasible, we enlarge N by adding more relay nodes in decreasing order of p * i , until the problem becomes feasible. Figure 5 shows in the red dashed curve the resulting minimum sum of powers obtained by setting P th = 1. Figure 6 shows the corresponding size of the set N . We observe in Figure 6 that the number of nodes with allocated power p * i exceeding P th = 1 (possibly including more relays to make the problem feasible) remains fairly constant as the size of the network n increases. This means that most of the power is allocated to a small subset of the nodes.
Removing the remaining nodes from the network and optimizing the power allocation again over the resulting simplified network results in the minimum total power plotted in Figure 5. This figure shows that even though the number of nodes in the simplified network remains approximately constant as n increases, the total power required to approximately achieve R 0 = 4 decreases with n. This is due to the fact that larger n allows to choose the best relays for the simplified network. There is some performance loss in terms of total power due to network simplification but this loss may be compensated by power savings arising from turning off some of the relays. While we have not modeled the power consumption of the clocks, CPU and other subsystems required to keep a relay active, in practice they may become comparable to the power consumed by radio transmissions. This makes network simplification very useful in practice.


Quantitative Comparison Between Crowd Models for Evacuation Planning and Evaluation

Crowd simulation is rapidly becoming a standard tool for evacuation planning and evaluation. However, the many crowd models in the literature are structurally different, and few have been rigorously calibrated against real-world egress data, especially in emergency situations. In this paper we describe a procedure to quantitatively compare different crowd models or between models and real-world data. We simulated three models: (1) the lattice gas model, (2) the social force model, and (3) the RVO2 model, and obtained the distributions of six observables: (1) evacuation time, (2) zoned evacuation time, (3) passage density, (4) total distance traveled, (5) inconvenience, and (6) flow rate. We then used the DISTATIS procedure to compute the compromise matrix of statistical distances between the three models. Projecting the three models onto the first two principal components of the compromise matrix, we find the lattice gas and RVO2 models are similar in terms of the evacuation time, passage density, and flow rates, whereas the social force and RVO2 models are similar in terms of the total distance traveled. Most importantly, we find that the zoned evacuation times of the three models to be very different from each other. Thus we propose to use this variable, if it can be measured, as the key test between different models, and also between models and the real world. Finally, we compared the model flow rates against the flow rate of an emergency evacuation during the May 2008 Sichuan earthquake, and found the social force model agrees best with this real data.

Introduction
Crowd simulation is an area that has been the subject of a significant amount of multidisciplinary work over the last few decades [1][2][3]. Its applications range from simulating crowds for movies [4,5] and games [6,7] to analyzing pedestrian behavior [8][9][10] and preparing for fire evacuations and similar emergencies [11][12][13]. The earliest attempts to simulate crowds generally adopted a macroscopic approach [14,15], where there is no explicit notion of an individual. Later, with increasing computational resources and with availability of observational data on an individual level [16][17][18], modelers were able to develop microscopic approaches [5,8] for application in areas where it was necessary to model and analyze individuals in the crowd. For example, in a simulation of evacuation, knowl-edge of the movement of crowds could reveal methods to improve crowd flow and evacuation speed.
One of the earliest and seminal works in individualbased motion planning was Craig Reynolds's model of coordinated animal motion such as bird flocks and fish schools [5]. Okazaki and Matsushita [19] assigned magnetic poles to goals, agents and obstacles to model movement. Subsequently, Helbing's Social Force Model [8] was developed, which is still one of the most popular models of movement in crowds. More recently, there have been several velocity-based approaches to motion planning, such as the synthetic vision based model [20] and the Reciprocal Velocity Obstacle Model [21][22][23].
The advocation of simulation-based analysis has become increasingly common over the last decade. Some well know applications include analysis of the yearly Mus-arXiv:1401.0366v2 [cs.MA] 10 Jan 2014 lim Hajj [24], or more recently the Love Parade disaster, Germany 2010 [25]. In the case of the latter, models and expertise had been used to guarantee the safety of the event, only for unforeseen circumstances to result in the deaths of 21 individuals. Clearly these real world examples emphasize the critical role that crowd modelling plays in safety preparation and planning, this in turn emphasizes the need for understanding the model dynamics, limitations and similarities. The extent to which these models are capable of accurately predicting the motion of a crowd is therefore critical for planning and safety.
Ideally these models should be validated against realworld data from all scenarios, for all cultures and for all varieties of crowd composition. Unfortunately real world data regarding egress or emergency situations is limited, often incomplete (the initial conditions are hard to know) and certainly not controlled. Even in the rare circumstances where data is available this often describes measurements and phenomena at the macro-scale, e.g., flow, density, average speed, etc. Because these models of individuals exhibit emergent behaviors at the mesoscopic and macroscopic scales, it is in general hard to tell whether the dynamics at the individual level are correct. The best researchers can realistically hope for is to collect microscopic data for a single scenario to calibrate the model. An alternative approach to the data-intensive one is to instead look quantitatively at the fundamental dynamics of the models and understand where these models differ, at the same time identifying common aspects of the models. This will aid in understanding fundamental aspects of all crowd models and offer insight into real-world dynamics of human crowds.
While much of the existing research does offer basic comparison of proposed models with existing models to demonstrate their usefulness, there is no existing quantitative comparison of the differences and similarities of these models to tell which is most accurate. The objective of this paper is to demonstrate a methodology for quantitative comparison of simulation models in a simple egress simulation. For this, we choose three popular models that are structurally very different. Since our focus is on the comparison methodology, we do not worry about specific versions of the models, in particular recent modifications and enhancements that are supposed to produce more realistic crowd behaviors.
The contribution of the paper is then three-fold: firstly, the analysis and comparison of these models provides interesting insight into the consequence of adopting each in particular forms of crowd or egress simulation. Secondly, the systematic approach we describe could be used in future to compare further models and develop a standard method of comparison for crowd simulation. Finally, the paper conclusion identifies a single measurable metric that is most effective in distinguishing the behaviour of the models.
The remainder of this paper is organised as follows. The models that are consider in the comparison are first described in Section 2; Section 3 describes the experiments that were conducted and the methodology for anal-ysis. Significant observations from the simulations are presented and analyzed in Section 4. Finally, Section 5 concludes the paper.

Models
In this section we describe the three individual (or microscopic) models that are compared in this paper. We use the general term agent to refer to individuals within each of the models. The lattice gas model is a probabilistic approach where the future location of an agent is probabilistically determined based on the current configuration of its neighborhood. This implies that the same initial configuration can produce different results during different runs. On the other hand, in the social force model and reciprocal velocity obstacle model, agents find their ways to their destinations via deterministic collision-avoiding calculations. As a result these models produce the same result on multiple runs for a particular initial configuration.

Lattice gas model
A Cellular Automation (CA) model is one in which space and time are discrete. Furthermore, the state space of a CA model is also discrete and finite. In each time step the values of all cells are updated synchronously based on the values of cells in their neighborhoods. Depending on the type of neighborhood (i.e., von Neumann, Moore), and the type of lattice (triangular, square, hexagonal, etc.), the exact number of cells in the neighborhood of a given cell can vary [26]. Lattice gas models are CA models that make use of a discretized version of the Boltzmann transport equation to model motion [27][28][29]. Advances in this modeling area include the Extended Floor Field Model [30] in which agents interact through virtual traces that act like the pheromones in chemotaxis, and the SWARM information model [31] which uses multiple floor fields to model transmission of knowledge between agents. In this study we use the model by Tajima and Nagatani [32]. In contrast to later CA models with additional features like bi-directional movement [33] and vision impairment during evacuation [29], this original model is simple, and adequate for our model comparison purpose.
A square grid is used, and each cell can be occupied by at most one person. The person performs a random walk biased towards the single exit in the room. For example, in Fig. 1(a), the probability that the person takes a step in the +y direction is P y = (1 − D)/3 + De y /(e x + e y ), while the probability that it takes a step in the −x direction is P −x = (1 − D)/3 + De x /(e x + e y ), since the intended direction e has positive projections along the +y and −x directions. Along the +x direction, for which e's projection is negative, the probability that the person takes a step in the +x direction is P x = (1 − D)/3. Here we see that the random walk probabilities contain an unbiased component (1 − D)/n, which is the same for all n permissible directions (in the Tajima-Nagatani model, the −y Intended direction

Intended direction
Intended direction direction is not permitted), as well as a biased component D = 0.7 that favours movement towards the exit.

Social force model
The social force model is one of the most popular models for motion planning in crowds [34][35][36]. This model is based on the idea that pedestrians move in response to fictitious attractive or repulsive social forces produced by obstacles and other pedestrians. Over the years, several extensions have been made to the social forces model like the modeling of grouping behavior [34]. In this paper, we use the Helbing-Molnár-Farkas-Vicsek (HMFV) social force model [8], which is tweaked from the original model introduced by Helbing and Molnar in 1995. In the original model proposed by Helbing and Molnar [8], three types of forces act on a given agent i with mass m i , instantaneous position r i (t), and instantaneous velocity v i (t). The first is a restoring force that steers the agent towards the desired velocity v 0 at a rate determined by the characteristic time τ i = 1. Here where V 0 is the preferred speed, e i is a vector that points towards the exit, and (1 − p) is the weight given to this desired velocity. With a weight of p, agent i also adapts to the average velocity v j i in its neighborhood. When p is small, agent i moves more along its intended direction e i (t), whereas if p is large, agent i tends to follow where its neighbors are going. We can therefore tune p from p ≈ 0 (self-directed normal egress) to p ≈ 1 (panic-driven herding during emergency evacuations). In this paper, we used p = 0.2 to simulate an emergency evacuation situation. The second is a repulsive force where that mimics the psychological tendency of agents i and j to move away from each other if they are too close.
Here R ij = R i + R j is the sum of radii of the two agents, d ij = |r i − r j | is the physical distance between the two agents, and n ij = (n 1 ij , n 2 ij ) = (r i − r j )/d ij is the unit vector pointing from agent j to agent i. Additionally, A = 2000 N and B = 0.08 m are the repulsion coefficient and the fall-off length of interacting agents respectively [37]. Helbing and Molnar also found it necessary to introduce two other terms in the interaction force when agents i and j are in contact with each other, i.e. d ij < R ij . This counteracting body compression term kη(R ij −d ij )n ij and sliding friction term κη(r ij − d ij )∆v t ji t ij are crucial for getting realistic behaviors of panicking crowds. Here t ij = (−n 2 ij , n 1 ij ) is the tangential direction and ∆v t ji = (v j − v i ) · t ij is the tangential velocity difference.
Finally, in the third force the first and second terms repels agent i from a wall that it is d iW away from, while the third term (which is negative) is introduced to mimic the observation that people move faster near walls when they are in crowded situations. Here, n iW is the normal vector of the wall, and t iW is the tangent vector of the wall. Also, k = 12, 000 kg/s 2 and κ = 24, 000 kg/ms are respectively the body force constant and the sliding friction force constant used.

Reciprocal velocity obstacles model
The third model we consider is the reciprocal velocity obstacles (RVO) model. In this model, the time to the next collision is calculated based on the relative velocities of N agents. Each agent then changes its velocity {v i } N i=1 to maximize this time to collision. By continuously updating the velocities in this manner, collisions are avoided. This algorithm was first proposed by Fiorini et al. [38], but was first used in multi agent systems in [39]. Since then there have been several modifications and improvements to the algorithm [10,22,40,41], although the underlying idea remained the same. In this paper, we use the RVO2 model introduced by Guy et al. [22], where the collision avoidance computation is based on computational geometry and linear programming.
Given a preferred velocity, RVO2 helps an agent find the velocity closest to the preferred velocity that will enable it to avoid collisions with all other agents. This is done by determining for each neighboring dynamic and static obstacle an Optimal Reciprocal Collision Avoidance line (ORCA line) . Each ORCA line determines the half plane in the velocity plane where the agent's velocity should lie to ensure that no collision occurs with that particular obstacle for the next τ seconds. τ is a parameter that specifies the number of seconds for which the chosen velocity should avoid a collision. From the set of half planes thus obtained, the optimal collision-avoiding velocity can be determined as the velocity within the permissible region for all the half planes that is closest to the preferred velocity. This can be determined efficiently using linear programming. Besides the speed and efficiency of the algorithm, another strong appeal of this model is that we need only set one parameter τ . For the experiments in this paper we use two values of τ ; 0.5 seconds or 10 simulation time steps for avoiding dynamic agents and 0.05 seconds or 1 simulation time step for avoiding static obstacles. The value of τ translates in practical terms to how early a walker tries to avoid a collision with an obstacle.

Methodology
In this section we explain the simulations that were carried out and the analysis done. Java 6 and the MASON simulation framework [42] were used for creating and running the simulations. The simulated scenario consists of agents evacuating from a simple rectangular room with a single exit. The dimensions of the room are shown in Fig. 2. Evacuation of the room was simulated with N = 50, 100, 150, 200, 300, 500 and 1000 agents whose locations were randomly distributed within the room. For each N , we ran 100 simulations for each model, and the positions and velocities of all agents at all time steps were collected. For meaningful comparison between the three models, we assumed that each agent is a perfect circle with radius R = 15 cm [43], that has a preferred speed of 1.3 ms −1 , and a maximum speed of 2.6 ms −1 [44]. For the social force model we also assumed that all agents have the same mass of 60 kg.
Where possible, we have used the same settings for all three models to ensure a fair comparison. The values used for the remaining parameters required for the models was given in Sec. 2 along with the respective models. However, to the best of our knowledge, only the social force model has been calibrated against real world data [45]. The other two models are not calibrated beyond choosing a reasonable and believable model and the values chosen are taken such that they produce empirically believable results.

Measurements
Depending on the context and motivation of the model there are several metrics that are generally used. For CA models that are generally used for studying macroscopic patterns, it is common to measure macro variables like the mean evacuation time [29] and flow rate [32]. For models like social force where the motivation is to study interactions at a more granular level, it is common to make empirical and qualitative observations like the lane effect [8]. Sometimes quantitative metrics like density-dependence of the flow or velocity are also measured [46]. From a graphics perspective, where performance is key, computation time [20], frame rate and run time per frame are generally measured against the number of walkers [23].
In this paper, we performed two main classes of measurements: time-based and distance-based measurements. For time-based measurements, we first measured the evacuation time distributions for different number of agents. To better understand the stages in the evacuation, we also divide the room into six different zones (see Fig. 3) to measure the evacuation time sub-distribution. The outer radii of Zones 1, 2, 3, 4, 5 and 6 are 5 m, 10 m, 15 m, 20 m, 25 m and 30 m respectively. We also measure the flow rate at the exit as a function of time for each model.
For distance-based measurements, we traced the trajectories of all agents to obtain the distribution of total distance travelled going from the initial position to the exit. Here, the total distance travelled D i by agent i is just the sum of its displacements Ti t=0 ||r i (t + 1) − r i (t)||, where T i is the evacuated time of agent i. We also divide the room into 100 × 100 cells, and count the number of agents passing through each cell as the simulation progresses. This is then visualized as a heat map. Finally, we measure the ratio of the total distance travelled D i by agent i to the minimum distance D min it would cover if it evacuated along a straight line (For lattice model, D min is the Manhattan distance between initial position of agent and exit). This in some sense quantifies the 'inconvenience', I i experienced by the agent i during the evacuation.

DISTATIS
After making six different measurements on three models, we discovered that when we look at the flow rates, the lattice gas model is more similar to the RVO2 model. However, if we look at the total distance travelled, the social force model is more similar to the RVO2 model. If we look at the zoned evacuation times, we find that the three models have very different distributions (see Sec. 4.1). Since we do not know a priori which measurements best discriminate the three models (or for that matter, best discriminate models from the real world), we want to be able to compare the three models quantitatively, incorporating information from all six measurements.
To do so, we compute the Jensen-Shannon divergence [47] D between the distributions f kµ (z) and f kν (z) for measurement k of models µ and ν. Here, z is a continuous variable like evacuation time, distance travelled, or inconvenience, is the Shannon information function, and If f kµ (z) and f kν (z) are highly similar, we will get D JS [f kµ , f kν ] ≈ 0, whereas if f kµ (z) and f kν (z) are very different, D JS [f kµ , f kν ] 0, i.e. the Jensen-Shannon divergence qualifies as a distance metric. In this way, we obtain six 3 × 3 distance matrices D k .
Then, we use the DISTATIS method [48] for analyzing multiple distance matrices. This is a generalization of the method of principal component analysis (PCA). STATIS is an acronym for the French expression 'Structuration des Tableauxà Trois Indices de la Statistique', which roughly means 'structuring three way statistical tables'. The difference between a straightforward application of the PCA method and the DISTATIS method is shown in Fig. 4. Instead of six different PCAs for the distance matrices obtained from the six different variables measured, and inevitably end up with conflicting conclusions on which models are more similar, we standardize the six distance matrices, and ask which variables are more similar to each other.
The idea behind DISTATIS is that the data points are similarly clustered in two variables, if their standardized distance matrices in these two variables are similar to each other. Therefore, if we compute the cross correlations between variables, we will discover which variables give more similar outcomes to which other variables through a PCA. Components of the first eigenvector obtained from this PCA tells us how important each variable is. By weighting the distance matrices of each variable by its component in the first eigenvector, we construct a compromise matrix whose matrix elements give us the most reliable similarity/dissimilarity between data points. A final PCA of this compromise matrix then gives the most reliable groups of data points. In Fig. 12 Fig. 4. In many classification problems, we can define different distance metrics for the same set of data points to measure how dissimilar they are to each other. In this figure we show three different distance matrices D1, D2, and D3. If we perform principal component analysis (PCA) on them separately (PCA1, PCA2, PCA3), we are not likely to agree on which data points are more similar to each other. In the DISTATIS method, we first standardize the distance matrices by centering and rescaling, to get S1, S2, and S3. If two data points are indeed more similar to each other than they are with the third, then we expect this similarity structure to be reflected in S1, S2, and S3, i.e. two standardized distance matrices would have very similar patterns of matrix elements. With this in mind, we reshape the matrices S1, S2, and S3 to make them into vectors of distances. We then compute the cross correlations between these vectors, before performing PCA on the correlation matrix C. This PCA tells us which distance matrices agree with each other more, and which less. The result is a consensus between different distance metrics, based on which a single PCA will discover the similarity structure between the data points.

Evacuation time
As Fig. 5 shows, the evacuation time distributions for all models have a uniform sub-distribution in the middle of the evacuation, which signals congestion at the exit, once we have N > 50 agents. However, in the RVO2 model this uniform sub-distribution is flanked by two small peaks. These two small peaks suggest a higher evacuation rate just before congestion sets in and just before congestion dissipates in the RVO2 model. This is not observed in real crowds, and hence is clearly an artifact of the model.
Comparing the three models, we find that the social force model has the longest evacuation times of up to five minutes. More importantly, careful inspection reveals fluid-like crowd dynamics around the exit in the RVO2 and lattice gas models. These can be seen in our videos http://www.youtube.com/watch?v= P64p3nlH_P4 and http://www.youtube.com/watch?v= vJ0Nzi5Bykw. In the RVO2 model, an obvious back flow can also be seen during the sudden gush of agents toward the exit. Due to the strong repulsion forces, we see a solid-like behavior in the social force model when the exit becomes congested (see our video http://www.youtube. com/watch?v=3Ot7m959_yo).
In addition, we divided the room into six zones, as shown in Fig. 3, based on the distance to the exit, to determine how strongly the crowds mix in the three models. The distributions of zoned evacuation times are shown in Fig. 6. These tell us that mixing is strong in the lattice gas and RVO2 models, but weak in the social force model. This mixing dynamics can be seen more clearly from our videos http://www. youtube.com/watch?v=qeoJotgEUxk (lattice gas), http: //www.youtube.com/watch?v=uZpd5LODYZs (RVO2), and http://www.youtube.com/watch?v=wEB6Ya0o_yw (social force). From the movies of the lattice gas and RVO2 simulations, we see that agents from the nearer zones get pushed to the side once there is congestion at the exit. In these two models, agents prefer to keep moving when it is possible for them to do so. Agents in the social force model behave differently: once the exit becomes congested, they will become nearly stationary and wait for their turn to go through the exit.

Passage density
The passage density heat map of agent locations averaged over 100 simulations provides information about the spatial-temporal trace of agents during their evacuation. This is crucial when analyzing the spatial structure of congestions. In Fig. 7(a), we see a channel leading straight through the exit. This channel acts as an attractor, because once an agent gets onto this, it will be forced statistically towards the exit (P t,y = 1.00).
Though not as pronounced as for the lattice gas model, two left-right symmetric channels also formed in the RVO2 model (Fig. 7(c)). These two channels point roughly from the centers of mass of the left and right halves of the room towards the exit. No such artificial channels can be seen in the passage density map of the social force model (Fig. 7(b)).
Comparing the three models, we also see that the RVO2 model (Fig. 7(c)) has the most compact passage density, whereas the social force model (Fig. 7(b)) has the least compact passage density. In particular, in the region just in front of the exit, the RVO2 model gives rise to very high passage density. This high passage density cannot be attributed to agents stepping over a cell once or twice, as in the lattice gas and social force models. In the RVO2 model, this high passage density is produced by repeated visits to the same cell by individual agents that are always moving.

Total distance traveled distribution and inconvenience
To measure the difficulty (inconvenience) for an agent to squeeze through the crowd during evacuation, we measured the total distance traveled and as well as the inconvenience (defined in Eq. 7). As we can see in Fig. 8(b) and Fig. 8(c), the distributions of total distance traveled are very similar for the social force and RVO2 models. The lattice gas model, on the other hand, has very different distributions of total distance traveled ( Fig. 8(a)). This is because space is continuous in the two former models, but discrete in the lattice gas model.
However, when we compare the distributions of inconveniences in Fig. 9, we find the lattice gas model is more qualitatively similar to the RVO2 model, though the typical inconveniences in the lattice gas model are larger than the those in the RVO2 model. For different N , the distribution peaks around the same inconvenience for both models. In contrast, the peak position of the inconvenience distribution of the social force model depends strongly on the number of agents in the simulation.

Evacuation rate
We also measured the flow rate of agents at the exit. As shown in Fig. 10(a), the flow rate of the lattice gas model saturates at around 11.8 persons per second due to the finite sizes of the agents and the exit. The social force model has the lowest flow rate of 2.2 people per second ( Fig. 10(b)), because agents slow down dramatically at the end as a result of strong repulsive forces between agents. From Fig. 10(c)), we see the same peaks before and after the congestion period for the RVO2 model, a consequence of the liquid-like dynamics of RVO2 crowds.
In an attempt to compare simulation results with real data we obtained a number of videos showing reallife evacuations [49,50] during the Sichuan Earthquake on May 12, 2008. We analysed two such videos, http: //www.youtube.com/watch?v=-y38ebiAnQw and http: //www.youtube.com/watch?v=e1yapP3z_L4. Because the geometry of the scene and the viewing angles are not known, we could only measure the flow rates on a secondby-second basis from these two videos (Fig. 11), for comparison against what we measured in our simulations. Unlike the computational study, where we could average over 100 simulations to get smoothly varying flow rates, the flow rates measured from the videos are noisy. However, sensible comparisons can still be made. In http://www. youtube.com/watch?v=-y38ebiAnQw, the crowd is shown to evacuate through the check-out counters of a Walmart supermarket. The flow rate was high, but there was no congestion, because of the wide exit. This flow rate cannot be compared against our simulations of evacuation through a narrow exit.
In http://www.youtube.com/watch?v=e1yapP3z_L4, students were evacuating their classroom through a narrow door. From the video, we see that the width of the door allows no more than three school children to exit simultaneously. Therefore, the exit dimensions seem to correspond to the exit dimensions use in our simulations. From Fig. 11, we see that the average flow rate in the congestion phase is about 3 persons per second. This is far below the congestion flow rate in the lattice gas model, and very close to the congestion flow rate of the social force model. Based on this observable, we find that the social force model agrees best with the real data. This small test also tells us that the three models can be distinguished through their congestion flow rates, if other observables could not be measured.

DISTATIS comparison
Finally, we make use of the DISTATIS method explained in Sec. 3.2 to compare the three models with all six different measurements. Fig. 12 shows the different measurements of the three models projected onto their first two principal components. The barycenters for the three models are also shown. From Fig. 12 we see that based on the total distance traveled (green), the social force and RVO2 models have similar distributions, which are both very different from that of the lattice gas model. On the other hand, the lattice gas and RVO2 models have similar distributions that are significantly different that of the social force model.
More importantly, we find the clustering of some observables: evacuation time, passage density, and flow rate all tells us that the lattice gas model is similar to the RVO2 model, and different from the social force model. We also find the inconvenience and zoned evacuation time provide some indication that the three models are qualitatively different. In fact, the DISTATIS method indicates that    the zoned evacuation time is one key metric that clearly distinguishes the three models. Therefore, any real world experiments should attempt to measure this. This observable maximally discriminates between the three models and therefore it is reasonable to expect it will also discriminate between real-world data and the models.

Conclusion
To conclude, we simulated emergency evacuation from a rectangular room with a single narrow exit using three models: (1) lattice gas, (2) social force, and (3) RVO2, with different starting number of agents in the room. From these simulations, we measured six observables: (1) evacuation time, (2) zoned evacuation time, (3) passage density, (4) total distance traveled, (5) inconvenience and (6) flow rate. Besides qualitative comparison of the various distributions obtained, we also compared the three models quantitatively using the DISTATIS method. Comparing the evacuation time, passage density, and flow rate distributions, we find that the lattice gas and RVO2 models are similar, both of which are very different from the social force model. On the other hand, if we compare the distribution of total distance traveled, the social force and RVO2 models are more similar to each other, and very different from the lattice gas model. We compared the simulated congestion flow rates of the three models to congestion flow rates of school children evacuating from their classroom during the May 2008 Sichuan Earthquake, and find based on this observable that the social force model agrees best with the egress data. Finally, from an analysis of the DISTATIS we have identified the zoned evacuation time as the one observable metric that can best discriminate between these models, and also between models and real-world data