Distributed Dynamic Quantile Solver with Plug-and-Play Operation

In this paper, we propose a continuous-time distributed algorithm for the dynamic quantile problem. The problem is to find the quantile of time-varying signals in a network of agents, each of which having the signal of its own. For example, this problem includes finding the median, maximum, or the second largest value of the signals. The proposed algorithm guarantees convergence from arbitrary initial conditions and does not use the decaying gains. Hence our algorithm is suitable for plug-and-play operation, where agents may freely join or leave the network during the operation. An application to a simplified electricity market problem is presented to show the effectiveness of the design.


I. INTRODUCTION
C ONSENSUS algorithms have been widely studied due to their distributed nature and simple operation [1]. In particular, estimating an average of given data, which is called average consensus [2], [3], is a well-studied problem for the multi-agent systems. Although the average is used to measure the central tendency of the data, it is not necessarily a robust measure since it can be easily biased by a few outliers.
An alternative statistical metric is a median which separates the higher half from the lower half of the data. Along with the median, maximum and minimum are other metrics measuring various statistical properties of the data. More generally, metrics mentioned so far can be represented by the quantile. The quantile is defined as a value that separates the higher percentages of the data from the rest, e.g., the median can be considered as the 0.5-quantile of the data.
There have been many studies to compute the aforementioned metrics in a distributed manner. However, most of the distributed algorithms for finding the maximum [4], [5], the minimum [6], [7], or the median [8] require a specific initialization process. Such an initialization process causes difficulties when considering a plug-and-play operation of agents in a network, that is, when some agent is added or removed during the operation of the algorithm. Specifically, if the network varies with time, such algorithms must be reinitialized for each change in the network, and this becomes a hindrance to the applications that require uninterrupted operations in spite of joining or leaving agents such as power network and V2G (vehicle-to-grid). The plug-and-play property for the multi-agent systems has attracted attention recently, and its necessities are well summarized in, e.g., [9], [10] (where the multi-agent systems with the plug-and-play property are called 'open' multi-agent systems).
Distributed algorithms computing the quantile without initialization process have been developed [11], [12]. However, most of the algorithms are based on decaying step sizes. The decaying step size not only slows down the convergence rate but also leads to problems for the plug-and-play operation. The reason is that whenever a new agent joins (or leaves the network), it takes longer to converge to the new solution due to the diminished step size compared to its initial value. Meanwhile, alternative solutions are proposed in [13]- [15] which solve the quantile problem in a distributed manner. However, subgradient descent based algorithm in [13] uses decaying step size, so it is hard to handle the time varying signals and plug-and-play operation. Besides, ADMM based algorithm in [13] does not use decaying step size, but it requires each agent to have additional memory space proportional to the number of its neighbors. Algorithms in [14], [15] VOLUME 4,2016 are limited to computing the median. More importantly, most of the algorithms mentioned so far only solved the consensus problem when the data are constant. It will be seen in the application of Section IV that the computation of quantiles for time-varying signals is sometimes necessary.
In this paper, we propose a distributed algorithm for the dynamic quantile problem, i.e., finding the quantile of timevarying signals. In particular, the proposed algorithm does not rely on initial conditions, which enables the plug-andplay operation. The proposed algorithm also does not use the decaying step size so that the convergence rate is linear, and the convergence time is explicitly guaranteed. The tradeoff is that the proposed algorithm only converges to the neighborhood of the desired solution. Nonetheless, we show that the approximation error can be made arbitrarily small by a suitable choice of parameters; moreover, the quantile can be obtained exactly in a finite time if all the signals are constant and the smallest unit of the constant signals is given.
This paper is organized as follows. Section II presents preliminaries on the dynamic quantile problem, and Section III provides our main results on the distributed algorithm. Then, Section IV presents an application to an electricity market problem. Finally, Section V concludes the paper.

A. NOTATION
The 2-norm of a vector x is defined by ∥x∥ := √ x ⊤ x and the induced 2-norm of a matrix A is written by ∥A∥. For a set M, the distance from a vector x to M is defined by ∥x∥ M := inf {∥y − x∥ : y ∈ M}. The absolute value of a scalar s is written by |s| and card(S) denotes the cardinality of a finite set S. For a real number s, the rounding up (resp. rounding down) is denoted by ⌈s⌉ (resp. ⌊s⌋). For scalars a 1 , . . . , a k ∈ R, we denote [a 1 · · · a k ] as a row vector. For matrices A 1 , . . . , A k , diag(A 1 , . . . , A k ) denotes the block diagonal matrix.
A coupling graph G := (N , E) is an undirected graph consisting of a node set N := {1, · · · , N } of nodes and an edge set E ⊆ N × N of ordered pairs of nodes, called edges, such that (i, j) ∈ E implies (j, i) ∈ E. In this paper, each agent is represented as a node of a graph and the information flow between two agents is represented as an edge, e.g., (i, j) ∈ E means that agent i gives its information to agent j. An adjacency matrix A := [a ij ] ∈ R N ×N is defined such that a ij = 1 if (j, i) ∈ E, and a ij = 0 otherwise. For a given node i ∈ N , N i denotes the set of all nodes j ∈ N such that (j, i) ∈ E. The Laplacian matrix L := [l ij ] ∈ R N ×N of a coupling graph G is defined as l ii = k̸ =i a ik for all i ∈ N and l ij = −a ij for all i, j ∈ N such that i ̸ = j.

II. PRELIMINARIES ON DYNAMIC QUANTILE PROBLEM
In this section, we define the p-quantile for time-varying signals and introduce the dynamic quantile problem. Then, we formulate the dynamic quantile problem as an optimization problem.
Consequently, the dynamic p-quantile problem is to find the value of the p-quantile of Z t at every time t. It is noted that if pN is not an integer, then ⌈pN ⌉ = ⌊pN + 1⌋. Hence, the p-quantile of Z t becomes a singleton.
The dynamic p-quantile problem can be formulated into an optimization problem. For a given 0 < p < 1, suppose that t is a fixed time and consider the following optimization problem where Example 2. For the case of Example 1, the function F p (s) is plotted in Fig. 1(a) with p = 0.3. It is seen that F p (s) is a convex and piecewise-linear function. Graph of the objective function J t (s) is shown in Fig. 1(b), where p = 0.3 and Z t = (1, 3, 6, 6, 8). It can be confirmed that the objective function is also a convex and piecewise-linear function. Moreover, the set of minimizers of the objective function is a singleton {3}, which is exactly the 0.3-quantile of Z t . ♢ Motivated by these observations, the following lemma shows that the solution of the optimization problem (1) at time t is the p-quantile of Z t . Lemma 1. Let M p,t be a set of minimizers of the optimization problem (1) at time t. Then, The proof is provided in the Appendix.

III. DISTRIBUTED DYNAMIC P -QUANTILE SOLVER
Although the optimization problem (1) solves the dynamic p-quantile problem as shown in Lemma 1, it is solved in a centralized fashion, i.e., signals of all agents must be gathered by some centralized unit to solve the optimization problem. On the contrary, distributed computation has many benefits (which are discussed in, e.g., [16]), and thus, we present a distributed algorithm, where each agent has its own signal z i (t) and communicates with its neighbors in the coupling graph G. The proposed distributed p-quantile solver is given byẋ for all i ∈ N , where α, γ > 0 are the design parameters and sgn(·) denotes the signum function such that sgn(s) = s/ |s| for s ̸ = 0 and sgn(0) = 0. It should be noted that each agent only exchanges the state x i and keeps the signal z i (t) privately, which can enhance privacy. The intuition behind the distributed algorithm (2) comes from the blended dynamics approach introduced in [17]- [19]. In this approach, we consider the heterogeneous multiagent systems whose dynamics are given bẏ where x i ∈ R is the state and f i represents the heterogeneous vector field of agent i. To analyze the behavior of the multiagent system, the blended dynamics of (3), which is the average of the vector fields of all agents, is introduced aṡ Then it has been shown that, if the blended dynamics satisfies suitable stability conditions, 1 then for any ϵ > 0, there is That is, the trajectory x i (t) is approximated by the solution s(t) of the blended dynamics.
The blended dynamics of (2) then becomeṡ which is the subgradient method for solving the p-quantile problem (1) when the signals z i (t)'s are constant. Hence, it 1 Details of the technical assumptions can be found in [17]. 2 It is the case when the zero vector is chosen as a subgradient of |x| at x = 0.
follows from (4) that x i (t) will converge to the neighborhood of the optimal solution M p,t as desired with sufficiently large γ. Consequently, even when z i (t) is no longer constant, we may expect x i (t) to be still close to M p,t if x i -dynamics (2) is relatively faster than the changing rate of z i (t). These observations are made rigorous under the following assumption. Assumption 1. The signals z i (t)'s are locally Lipschitz and there exists B ≥ 0 such that the signals satisfy 3 with Assumption 1, we obtain the following result.
Theorem 1. Suppose that the coupling graph G of all agents is connected and that Assumption 1 holds. Let the design parameters α and γ be such that Then, there exists a unique solution of the p-quantile solver (2) from any initial condition x i (0) ∈ R, and the solution satisfies where λ 2 is the second-smallest eigenvalue of the Laplacian matrix. ♢ Remark 1. The second-smallest eigenvalue λ 2 (i.e., the algebraic connectivity) measures the connectedness of the graph and depends both on the topology of the graph and the number of the agents N . In general, λ 2 is lower bounded by 4/(N 2 − N ) for all connected, unit weight, coupling graphs G [20]. Thus, if the upper bound of the number of the agents N throughout the operation is known, then we can design the parameters α and γ such that the steady-state error bound in (6) remains constant even under the plug-and-play operation, where some of the agents join or leave the network. Indeed, for any constant b > 0, let so that (5) holds for all 1 ≤ N ≤N , then the steady-state error bound in (6) becomes ♢ Remark 2. From Theorem 1, it follows that the trajectory of each agent x i (t) converges to the neighborhood of the optimal solution M p,t regardless of its initial condition x i (0).
From this, we can conclude that the proposed algorithm is suitable for the plug-and-play operation, in the sense that, when some of the agents join or leave the network so that the number of agent N and the network topology are changed (but still remains connected), Theorem 1 applies to the new situation and the proposed algorithm (2) guarantees the estimate converges to the neighborhood of the new optimal solution M p,t . We demonstrate this in Section IV. ♢ Proof. For the stacked signal z(t) := [z 1 (t) · · · z N (t)] ⊤ ∈ R N and the stacked state x := [x 1 · · · x N ] ⊤ ∈ R N , the overall system (2) can be written aṡ where L is the Laplacian matrix of the graph G, SGN(x) := sgn(x 1 ) · · · sgn(x N ) ⊤ , and 1 N : Since (7) is a discontinuous dynamical system, the existence of the solution is considered in the sense of Filippov [21]. Let Then, H(t, ·) is a convex function for fixed t, and the right-hand side of (7) is a negative subgradient of H(t, x) with respect to x at (t, x). Thus, the existence and uniqueness of the Filippov solution of (2) are guaranteed by [21, Proposition S2, Proposition S3, and Proposition 9]. Under a connected graph G, the Laplacian matrix L is symmetric and has N − 1 real positive eigenvalues and a simple zero eigenvalue. It is supposed, without loss of generality, that the eigenvalues λ i for all i ∈ N of the Laplacian matrix satisfy 0 = λ 1 < λ 2 ≤ · · · ≤ λ N , where λ 2 is called the algebraic connectivity. By the symmetry of the Laplacian matrix, there exists an orthogonal matrix U ∈ R N ×N such that where Λ + := diag(λ 2 , · · · , λ N ) ∈ R (N −1)×(N −1) and R ∈ R N ×(N −1) such that R ⊤ R = I, R ⊤ 1 N = 0, and 1 N 1 N 1 ⊤ N + RR ⊤ = I. Then, ∥R∥ = ∥R ⊤ ∥ = 1 and L = RΛ + R ⊤ . Let r i ∈ R (N −1) be the i-th column of R ⊤ , i.e., R ⊤ =: r 1 · · · r N . Then ∥r i ∥ = (N − 1)/N ≤ 1 for all i ∈ N . Now, we prove the ultimate boundedness (6). Consider a coordinate transformation in whichx = (1/N ) i∈N x i ∈ R is an average of the state and ξ = i∈N r i x i ∈ R (N −1) so that x = 1 Nx + Rξ, i.e., Then, (7) is transformed intȱ First, we prove ξ(t) is ultimately bounded. Let W (ξ) := ||ξ||. Then, we have the following two cases: (ii) If ξ(t) = 0, then W (ξ(t)) = 0, and In either cases, we obtain
By using the similar argument as in the first case, we obtain x i (t) < z * ⌈pN ⌉ (t) for all i ∈ N , card({i ∈ N : x i (t) < z i (t)}) ≥ ⌊(1 − p)N + 1⌋, and for almost all t, whereδ := ⌈pN ⌉ − pN ∈ [0, 1). Thus, for almost all t ≥ T 1 (ϵ), if either Case 1 or Case 2 holds, i.e., if V (t) ≥ (α √ N − 1)/(2γλ 2 ) + ϵ, then we have where β * := (α/N )(1 − max{δ,δ}) − B > 0. The strict positivity of β * follows from the design of α in (5). By the comparison lemma, it is seen that the upper bound of V (t) decreases with constant speed β * when t ≥ T 1 (ϵ) and V (t) is larger than (α √ N − 1)/(2γλ 2 ) + ϵ. Let where Then, T 2 (ϵ) is non-negative and decreasing function in ϵ. Hence, we obtain for all t ∈ [T 1 (ϵ), T 2 (ϵ)] and for all t ≥ T 2 (ϵ). The guaranteed upper bounds (14) and (15) of V (t) are shown in Fig. 2. Hence it follows that ∥x(t)∥ Mp,t is ultimately bounded. Finally, we prove (6). From (8), we have Let where Therefore, Remark 3. Instead of the asymptotic result (6), the proof of Theorem 1 also yields the following non-asymptotic statement: for any error bound ϵ > 0, it holds that Here, T • (ϵ) represents the guaranteed time for the distance from x i (t) to the p-quantile M p,t to be less than or equal to the steady-state error bound with ϵ-margin, which is computed from (11), (13), and (17). It depends not only on ϵ, but also on the initial conditions ∥x(0)∥ Mp,0 and ∥ξ(0)∥, because T • (ϵ) is defined by T 2 (·) which depends on the initial conditions. By (11), it is clear that T • (ϵ) tends to infinity as ϵ tends to zero. The convergence rate of x i is linear by (16) because the convergence rate ofx is linear (as seen by (14) and in Figure 2) while the convergence rate of ξ is exponential (as seen by (10)). (It is also noted from (11) that T 1 (ϵ) in Figure 2 can be made arbitrarily small by choosing γ sufficiently large.) ♢ Remark 4. The proposed algorithm has the freedom of choosing two design parameters α and γ. The relation among the number of agent N , the size of the error in (6), the convergence time T • (ϵ) with ϵ-error in (18), and two design parameters α and γ, can be also seen from the proof of Theorem 1. Indeed, if α is chosen as the order of α = O(N ) VOLUME 4, 2016 so that the assumption of Theorem 1 is met, and if γ is chosen as O(N 3.5 ), then the steady-state error in (6) and the time T • (ϵ) in (18) are both O(1). ♢ By Theorem 1, the proposed algorithm (2) approximates the p-quantile M p,t at time t, which is an interval in general. However, if pN is not an integer, the p-quantile is a singleton as mentioned in Section II.
It is often useful to estimate the value of the specific rank rather than the interval. The following corollary shows that the proposed p-quantile solver (2) can be used to estimate the k-th smallest signal z * k (t) if the value of p is appropriately chosen utilizing the information of k and N . Corollary 1. Assume (k − 1)/N < p < k/N for some k ∈ N , G is connected, and Assumption 1 holds. If the design parameters α and γ satisfy (5), then there exists a unique solution of the p-quantile solver (2) from any initial condition x i (0) ∈ R, and the solution satisfies ♢ Proof. By the assumption, we have k − 1 < pN < k. Then, the set of minimizers of the optimization problem (1) at time t is given by Hence, it follows from Theorem 1 that the solution exists uniquely and satisfies lim sup Although the steady-state error bound can be made arbitrarily small by adjusting the parameter γ, the solution of the proposed algorithm (2) only converges to the neighborhood of the p-quantile M p,t . However, the algorithm can compute the exact p-quantile in a finite time using the round-off operation if all the signals are constant (B = 0) in an evenly spaced grid ρZ := {ρm ∈ R : m is an integer} for some ρ > 0.

IV. APPLICATION TO ELECTRICITY MARKET PROBLEM
To appreciate the proposed algorithm, let us consider a simplified electricity market problem for a smart grid with mobile power sources, such as V2G (vehicle-to-grid). Suppose that many mobile power sources (e.g., electric vehicles) are connected to a power grid, and each source can sell and provide an amount of power generation, G, to the grid if a condition is met. The condition is that the sale price p i (t) of each source node i (which is decided by the source node itself) is less than or equal to the real-time market price. The price of each source may vary with time, as in Fig. 3, because the owner of power source may want to make more profit at peak times of electricity consumption. The power grid also includes the nodes that consume electricity. Let us suppose that each node can be a consumer as well as an electricity source. (This is for simplicity of simulation, but they can, of course, be different nodes, or their role can be switched depending on time. An electric vehicle is an example of such a node because they often sell the power that is stored overnight.) As a consumer, their power demand is denoted by d i (t) for node i, which is timevarying and decided real-time by the node.
The market for the power grid buys electricity from the source nodes of the cheapest price until the total demand is met. For this, the market controls the market price. For example, consider the case of D(t) = 200 and G = 100, for which two source nodes are to be identified. For this, the market price begins with the cheapest one from the sale prices {p i (t)} to select the first suitable source node, and then, increases the market price to the second cheapest one to identify the second source node. With the variation of D(t), the market price should vary in order for the market to buy from additional nodes or to drop the nodes, or even to change the nodes if the node for the cheapest sale price is changed. Since the node is mobile, some node may also leave the network or new nodes may join the network during the operation of the market for the power grid. The market price should be determined by considering all these aspects.
The problem is to decide the market price in a distributed manner for automatic selection of the source nodes. Our solution begins with a distributed estimator of total demand D(t): is the time-varying index set of neighboring nodes (which is available information to node i). The quantitŷ N (t) is an estimation of the number of participating nodes to the grid which is time-varying and can be estimated in a distributed manner by various methods such as [23]- [25]. Among them, we have used [24] in the simulation, which is suitable for our plug-and-play operation. As time goes by, each state x D i (t) approximately converges to the total demand D(t) (in particular, ∀ϵ > 0, ∃γ D such that lim sup t→∞ |x D i (t) − D(t)| ≤ ϵ for all participating nodes; for details, see [26]). Fig. 4 shows the states x D i (t) for all nodes, each of which approximately estimates the total demand.
Finally, the market price is distributedly determined by the proposed p-quantile solver: in which the quantity p of (2) is chosen as ( x D i /G − 0.5)/N (t) in the spirit of Corollary 1 with k = x D i /G ; for example, if x D i /G = 2 then x P i (t) approximately converges to the second cheapest price among {p i (t)}. Since each node runs the above algorithm, the market price becomes available to every sources, and so, whenever their sale price is less than or equal to the market price (i.e., p i (t) < x P i (t) + ϵ with small margin ϵ > 0 that is introduced because the proposed algorithm only approximates the true market price), the sources node provides its power G to the grid.
Simulation result is shown in Fig. 5, for which we used G = 100, α D = 160, α P = 500, γ D = 10000, γ P = 3000, ϵ = 0.5. It is clearly seen that, when D(t) is less than or equal to 100, i.e., ⌈D(t)/G⌉ = 1, the cheapest sale price is chosen as the market price (for the time from 0 to 6), and when 100 < D(t) ≤ 200, i.e., ⌈D(t)/G⌉ = 2, the second cheapest is chosen (for the time from 6 to 10). In the simulation, we used ode45 (Runge-Kutta method) for our algorithm, which is a built-in ODE solver in MATLAB. VOLUME 4, 2016

V. CONCLUSION
In this paper, we proposed a continuous-time distributed algorithm for solving the dynamic quantile problem. Finding the p-quantile is first interpreted as an optimization problem. Then the blended dynamics approach is used to obtain a distributed algorithm that converges to the neighborhood of the desired solution. We also characterized the relationship among the coupling gain, the steady-state error bound, and the convergence time. A main benefit of the proposed design is that it converges to the neighborhood of the p-quantile solution regardless of any initial condition and uses a fixed gain which allows for plug-and-play operation. The proposed algorithm is applied to the electricity market problem, and the simulation results verify the proposed algorithm.