Distributed optimization algorithms for game of power generation in smart grid

. In this paper, we consider a problem of ﬁnding optimal power generation levels for electricity users in Smart Grid (SG) with the purpose of maximizing each user’s beneﬁt selﬁshly. As the starting point, we ﬁrst develop a generalized model based on the framework of IEEE 118 bus system, then we formulate the problem as an aggregative game, where its Nash Equilibrium (NE) is considered as the collection of optimal levels of generated powers. This paper proposes three distributed optimization strategies in forms of singularly perturbed systems to tackle the problem under limited control authority con- cern, with rigorous analyses provided by game theory, graph theory, control theory, and convex optimization. Our analysis shows that without constraints in power generation, the ﬁrst strategy provably exponentially converges to the NE from any initializations. Moreover, under the constraint consideration, we achieve locally exponential convergence result via the other proposed algorithms, one of them is more generalized. Numerical simulations in the IEEE 118 bus system are carried out to verify the correctness of the proposed algorithms.


1.
Introduction. Recent decades have witnessed an explosive growth in electricity demand under limited energy sources causing a significant environmental and economic concerns for energy management system (EMS) in power grids. In order for utility companies to enhance energy management efficiency, there are a number of strategies have been proposed. These strategies include game-theoretic power consumption scheduling program and load billing [17], [5], advanced direct load program [22] and interruptible load schedule [21], to name only a few.
We study here a problem of controlling the power generation at electricity users (consumers) side, in which each user competes others to seek its optimal generation level in order to selfishly maximize its own benefit or minimize cost function while respecting distributed constraints as follows where x i is the produced power constrained in the set A i known only by user i, α presents global parameters or information of the grid,x ∶= ∑ i∈P x i and P is the set of users with its cardinality P . The problem (1) refers to an aggregative game, a special class of non-cooperative game, in which the individual cost function of each player (or player in the context of game) depends on its action and the sum of all players' actions, and the considered set of optimal levels {x * 1 , ..., x * P } constitutes a Nash equilibrium (NE) of the game [10]. As the power network's size increases, solving the optimization task (1) by centralized strategies might be impractical; hence there is a need for distributed algorithms. This motivates us to design distributed strategies exponentially converging to the optimal solution (1).
Our work finds relevance in the larger context of energy networked problems [17], [5], game theoretical works [10]- [31], limited control authority (LCA) [12], [13] and singularly perturbed system (SPS) [11]. The reference [17] studies a selfish optimization of each user's power cost, and proves that it is equivalent to the minimization of whole grid power consumption. To seek collective optimal solutions forming a NE in a game in a distributed fashion, each user executes interior point method. It, however, is assumed that each user directly knows detailed scheduled power consumption of others by exchanging information via a fully connected communication network. For tackling this drawback, the work [5] formulates the optimization problem in [17] as an aggregative game with a generalized billing scheme, then solves it by designing a distributed asynchronous gossip-based algorithm employed from [4]. This work only requires the average power consumption of the whole grids to be estimated by each user through communicating with its immediate neighbors. The authors in [5] and [14] share the similar idea that only one player is supposed to be wake up at each tick of a virtual clock to randomly communicate with one neighbor for estimating the time-varying sum of all players' decisions while forcing its action toward NE of aggregative games under constraints. The work [25] further extends [14] to more general games by excluding each agent's decision from its estimation vector on players' actions and taking constant step size into account. Another discrete-time algorithm, alternating direction method of multiplier (ADMM), is recently considered as an effective tool for dealing with distributed seeking NE problem [26]. Compared with a gossip-based method in [25], the designed inexact-ADMM in [26] performs faster convergence rate by defining an extra penalty term as an error between the estimates of each player's neighbors and using synchronous updates on players' actions and their estimate vectors. In a recent work [28], the authors propose linearized ADMM-like NE seeking approach using asynchronous updates and fixed step size on gradient descent while ensuring a fast convergence rate as synchronous updates.
None of the approaches mentioned above, however, studies the seeking NE of aggregative games under the LCA. For a network of n agents, in which each one governed by a distributed dynamic asẋ i = d i , where d i is the driving command, for i ∈ {1, 2, ..., n}, the LCA is naturally associated with physical limitations on agents' actuations [12]. Specifically, in dynamic average consensus (DAC) problem, when each agent has an authority to access time-varying input signal u i (t) ∶ [0, ∞) → R, and attempts to track the average of the network's inputs: x i → 1 n ∑ n j=1 u j (t), the LCA requires different slow desired convergence rates of agents towards the average if f i is implemented in real time based on a consensus algorithm. This requirement yields, for instance, is an initial value, meaning that the tracking rate needs to be less than or equal γ i > 0. The work [18] addresses the LCA when designing algorithm to deal with the problem (1). It, however, assumes the identity of convergence rate, i.e., γ 1 = γ 2 = ... = γ n , which might be impractical in networked systems. In contrast with the previous work [18], this work considers different convergence rates toward dynamic agreement under the LCA, hereby seeking the NE at each player's own pace. Our work herein is also related to SPS comprised of slow and fast motions, which often occur when we mathematically model some natural processes, for example, biochemistry, fluid dynamics, nonlinear kinetics, and so on [27]. We design algorithms in the form of SPS to seek NE of the game (1). Moreover, different from the aforementioned discrete-time algorithms [17], [5], [14]- [28], we develop distributed continuous-time algorithm to seek NE of an aggregative game under constraints. One of proposed methods in our work herein is similar to networked communication protocol of [28] and [31] where all players only agree to exchange message with their immediate neighbors rather than directly observe the remaining players' states as be assumed in [17]. This leads us to define an aggregate vector associated with each player to keep the time-varying estimation on other players' information.
Consequently, as the first contribution of this paper, we reformulate the selfish optimization problem of power generation in SG as a game. The physical layer is based on the infrastructure of IEEE 118 bus system, whereas distributed communication in the cyber layer is assumed possible due to smart meters at all units in the network, which helps update global information by distributed fashion. Remarkably, the model in this paper is a generalized one over of energy network model as we already mentioned in [18] and [19]. As our second contribution, we study the optimal energy generation from the view of game theory. It is proved in this paper that different desired rates toward dynamic consensus value among users do not effect to the NE. In comparison with other works in the literature [17], [5], [14]- [31], our work strongly touches on the LCA issue when designing distributed algorithms. To the best of our knowledge, the seeking NE issue under the concern of LCA has not been studied before. Our third contribution refers to the generalized algorithm in Section 6. Here, since each player can indirectly estimate other players' information by only communicating with its neighbors, this method is able to be applied for seeking NE of more general games.
In Section 2, some preliminary results in the SPS, limiting equation theorem and convex optimization are provided. We formulate the problem as a game in Section 3. Section 4, 5 and 6 are dedicated to present our distributed strategies. Numerical simulation are illustrated in Section 7. We conclude this paper in Section 8, where we discuss future research directions. All proofs are given in Appendix, together with a basic graph theory.
We denote 1 n ∈ R n as a column vector with all entries be 1, and 1 n×n ∈ R n×n be a matrix with all elements be 1. Matrix and vectors are in bold faces. The cardinality of set D is D . Let (D, d) be a metric space. A close δ-neighborhood of a set B, which is the set of all nonempty compact subsets of D, is defined by is the distance from a point x ∈ R n to a nonempty compact set B. The relative interior of the set D, denoted relintD, is defined as where aff (D) is the affine hull of a finite point set D in Euclidean space R n . The subscript r and q refer to reduced states and quasi-steady states. The terms G.A.S., G.E.S. and E.S. refer to globally asymptotically stable, globally exponentially stable and locally exponential stable, resp.

2.2.
Limiting equation theorem. We consider the nonlinear non-autonomous systemẋ = f (x, t), where f (x, t) = g(x) + d(x, t) with d(x, t) is a perturbation term to be guaranteed as lim t→+∞ d(x, t) = 0, and this system satisfies one of the following assumptions: (S1) f (x(t), t) is bounded for any bounded x.
If the trajectories of the unforced systemẋ = g(x) asymptotically approach the limiting set S and S is bounded, then this set is also asymptotically approached by all bounded solutions of the systemẋ = f (x, t) [2].

2.3.
Convex optimization with distributed constraints. We consider the problem of solving convex optimization associated with a network of agents in term of distributed inequality constraints.
where H ∶ R n ↦ R, h i ∶ R ↦ R are continuous differential and convex functions, m is a number of constraints. Let D be the intersection of the domains of all function h i . Regarding to Slater's theorem, if Slater's condition holds, i.e., there exists a strictly feasible point x ∈ relintD such that h i (x) < 0, i ∈ {1, 2, ..., m}, then strong duality holds. The Slater condition can be refined if some s of the functions h i (x) are affine and there exists x ∈ relintD with h i (x) ≤ 0, i = 1, 2, ..., s and h i (x) < 0, i = s+1, 2, ..., m. The strong duality implies that x * is the solution of (3) if and only if there exists λ * = [λ * 1 , λ * 2 , ..., λ * m ] ⊺ ∈ R m , together with x * , satisfies the Karush-Kuhn-Tucker (KKT) condition [3] : where λ * i h i (x * ) = 0 refers to complementary slackness condition. Specifically, we express strict complementary slackness condition if the following assumption is satisfied: 3. Game formulation. We model a network as a combination of n buses, where each bus belongs to one of four sets: a power load unit set L = {L 1 , L 2 , ..., L L } consisting of only electricity loads, a power generation unit set G = {G 1 , G 2 , ..., G G } comprising of only generators, a power load-generation unit set P = {P 1 , P 2 , ..., P P } presented for a set of electric users equipped with generators and considered as a set of players in the game, and finally a distribution transformers unit denoted as T = {T 1 , T 2 , ..., T T }. Let V = P ∪G∪L∪T , then n = V . Each unit is equipped with a smart meter integrated with a CPU and peripherals, thus each meter directly accesses to real-time data of its consumption or generation energy, and exchanges information with its neighboring meters via data links. The proposed distributed algorithms implemented in CPUs enable each unit to compute a consensus information in real-time, then based on that, each player calculates an optimal set point to operate its generator in order to get benefit selfishly in the game. The network of the smart meters, e.g. Local Area Network with Ethernet technology, can be modeled as a connected and undirected graph; hence, graph theory can be applied to analyze the interactions between them (see Appendix A). Remarkably, the model developed by graph theory concept can treat the transformers' losses and does not require any additional centralized communication channel to broadcast the information of the total consumed power of the grid to all users as in our previous works [18] and [19]. Without loss of generality, we suppose that the sampling time is slotted into hour-long intervals. At each slot, we respectively denote c L i , c G i , c P i , c T i as the consumed power of load i ∈ L, the produced power of generator i ∈ G, the consumed power of player i ∈ P , and the power loss at transformer i ∈ T . These values are positive constants and considered as uncontrollable variables, and also assumed to be exactly measured by smart meters at each slot of time. For the sake of analysis, we consider c i as one of c L i , −c G i , c P i or c T i if i ∈ L, G, P, T correspondingly. When the player i consumes a power, it needs to pay the charge. Thus, in order to curtail that monetary power expense, instead of only purchasing power, each player i needs to produce a power denoted as x i , selling it to the grid and getting some benefits. However, it comes at a cost that each player has to pay for power generation (e.g., fuel cost), yielding trade-offs between electricity bill payment and operational cost. Thus, each player is at most profit if it can seek its optimal power generation and operate at that level defined as x * i . From a game theory perspective, these optimal values of all players would form a NE, which is the point that none of the players (or decision makers) gets profit by unilaterally moving their actions given the actions of the others in the game. It means that if we denote x = [x i , x −i ] ⊺ ∈ R P as a strategy vector in which player i plays an action x i ∈ R and the remaining players play an action profile is an associated cost function of player i and α is a parameter. Specifically, in our game formulation, the NE is the solution of the following optimization problem The other term J i (x, α) refers to power price unit corresponding to billing scheme as in [17], [18] and [16], and is of form [19]). To solve (5), the global information ∑ i∈V c i −x needs to be precisely estimated at each player. We then trivially obtain Then under the assumption on convexity of the set A i , the game possess at least one NE [24].
Furthermore, since the Hessian matrix B = ∂ 2 fi(xi,x,α) ∂xi∂xj ∈ R P × P is symmetric matrix and ∂ 2 fi(xi,x,α) , one can derive that B + B ⊺ ≻ 0, leading to the unique existence of NE of the game [24]. In the following sections, distributed algorithms are presented to seek the unique NE, in which the sum 1 n (x − ∑ i∈V c i ) is estimated in distributed fast dynamics, and the NE is simultaneously calculated in slow dynamics. We address the problem that each player is adopted to the LCA, and thus it cannot track the dynamic consensus 1 n (x − ∑ i∈V c i ) at high rate but at its own slow feasible rate.

4.
A distributed algorithm for the game without constraints in power generation. In this section, we design distributed algorithms on the basis of the SPS in Section 2.1 and DAC protocol from [13]. In contrast with static average consensus (SAC), in which each agent reaches the mean of constant inputs, the DAC is applied when each agent tracks the mean of time-variant inputs. For the game without constraints in the power generation, we propose the Algorithm 1.
In the Algorithm 1, l ij is an entry of Laplacian matrix L (see Appendix A) and k > 0 is a constant. The action x i and its derivativeẋ i are considered as slowly time-variant inputs for the fast dynamics. When i is 0, the fast dynamics behave similarly as SAC in terms of the convergence of z i and w i toward agreement states z q i = −(1 n) (x − ∑ i∈V c i ) and w q i = −(1 n) ∑ i∈Pẋi , which are referred to quasisteady states in the SPS. In practice, the exchange of information in network under the LCA requires a certain amount of time, leading that each player can only track the agreement at its own positive pre-determined rate. Since the designed slow dynamics rely on the agreement value and gradient-based method to seek the NE, the rate of convergence towards NE is bounded above by the pre-determined rate. Associated with each player i, we define y i as its tracking value of the agreement

Algorithm 1 Without power generation limits
at a distributed rate γ i under the LCA. One can achieve a reduced equation for y i by substituting the quasi-steady states z q i and w q i into equation (8) Then, according to Theorem 11.2 [11], we have where κ i (0) = y i (0) − 1 n (x(0) − ∑ i∈V c i ). This motivates us to define y i as a value tracking 1 n (x − ∑ i∈V c i ) and to utilize the function We shift the system (7)-(8) to equilibrium points (z * , µ * , w * , ν * , y * ) corresponding to the NE x * by considering variable changes asẑ which can be verified to follow the condition (B1), andĉ i = 0 as c i is constant at each slot of time. To derive boundary-layer model, we rely on the fact that the fast dynamics (11) globally exponentially converge toẑ q i = − 1 n 1 ⊺ P x and globally converge toμ q i [8]. These states are isolated roots of (11a) when i = 0, satisfying the condition (B2). From (11a), by definingξ z we arrive at a boundary-layer model as follows and also obtain an analogous boundary-layer model for (11b)-(11d). A system constructed from (13) for ∀i ∈ P, j ∈ V is G.E.S and uniformly in (t,x), which is originally analyzed in [8] and thus omitted. Thus, we can obtain the following result  (11) is G.E.S., uniformly in (t,x).
We then study the convergence and stability property of reduced model for the transformed system (11)-(12) around the NE.  There exists * > 0 such that for all 0 < i < * , i ∈ V , the origin of the transformed system (11)-(12) is G.E.S. Therefore, if each player unit i ∈ P adopts the proposed distributed algorithm (7a)-(7b) and (8), whereas others follow (7c)-(7d), then its action eventually reaches the optimal power generation level in a distributed fashion from any initialization.
Proof. See Appendix B.

5.
A distributed algorithm for the game under constraints in power generation.
Lemma 5.1. The non-cooperative problem (5) is equivalent to the convex optimization problem with distributed inequality constraints. Here, , then H(x, α) is strictly convex and its optimizer x * is unique. Moreover, under the Slater's condition in Section 2.3, there exists λ * so that (x * , λ * ) satisfies the KKT condition (4), and the existence is unique according to Lemma 1 [18].
It is worth distinguishing the optimization problem (14) from distributed energy resources (DERs) [1] or economic dispatch problem (EDP) [6], in which all the agents in an interconnected networked system work cooperatively to minimize the global cost ∑ i∈P f i (x) under constraints. Our previous work in [18] was developed under the assumption that the identity of the feasible convergence rates of all players towards the agreement estimation are same, which yields a major restriction in practical issues. Therefore, in this work, we generalize the previous work by consider different convergence rates among them. For the game under constraints in power generation, we propose the Algorithm 2.

Algorithm 2 With respecting power generation limits
Fast Dynamic: (7) Slow Dynamic: for i ∈ Ṗ We will show that the trajectories of the dynamical system in the Algorithm 2 converge to the (y * , λ * , x * ), where x * is the unique NE of the game. For the sake of analysis, we denote the set composed of ζ elements, then ς + ζ ≤ P , and D A l ∩ D A u = ∅ because generator cannot operate at the minimum and maximum value at the same time. We also define , and assume strictly complementary slackness condition is satisfied, i.e., x * i > x l i for i ∈ D I l and x * i < x u i for i ∈ D I u , which leads to P = D A ∪ D I and D A ∩ D I = ∅. By following an analogous analysis as presented in Section 4, we obtain a transformed system (15) around (y * , λ * , x * ) as followṡ Moreover, the above transformed system in slow dynamic can be verified to satisfy condition (B1) and a reduced model is obtained aṡ ..,f ⊺ P (σ P , t)] ⊺ ∈ R 3 P . By following the analysis in Lemma 4.2, one can obtain a reduced model of the system in the Algorithm 2 aṡ wheref y (ŷ r ,σ r , t) = −kBŷ r + k n 1 P × P (λ r l −λ r u ) + kDdiag{δ r i (0)}r(t),λ r l = [λ r l1 , ...,λ r l P ] ⊺ ∈ R P ,λ r u = [λ r u1 , ...,λ r u P ] ⊺ ∈ R P . For the sake of analysis, we separate the right-hand side of two subsystems (18) asf y (ŷ r ,σ r , t) = f y (ŷ r ,σ r ) + are satisfied. There exists * > 0 such that 0 < i < * , ∀i ∈ V , if every unit i ∈ V adopts the proposed Algorithm 2 and the initialization of λ is positive, then each player's action eventually reaches its optimal power generation level locally in a distributed fashion.
Proof. See Appendix B.

6.
A generalized distributed algorithm for the game via leader-following consensus approach. The algorithm in the Section 5 is specifically designed for the aggregative game since only dynamic average of information is required to be estimated by all players. This restriction motivates our more general design, which can be applied for more general games, and the aforementioned game, therefore, is only one special case. In general problems, f i (x) can be a quadratic function [29], an aggregate of logarithmic functions in congestion game [30], or specific potential functions in potential game [15], etc. To tackle these problems, the idea we consider here is that each player can solve more general games if it knows other players' actions. Related to this topic, in [25], [28] and [31], each player utilizes an aggregate vector to store temporary estimation on others' information which is computed through distributed methods. But, different from these works, we consider the LCA that enables each player to track the information of others at its own pace.
We apply a leader-follower consensus approach in [23] to our design. Each i ∈ V is considered as a virtual leader broadcasting a reference information denoted as That reference cannot be observed directly by nonneighboring units but estimated via message exchanged among neighbors in the communication graph. To this end, we define a local vector associated with player i as y i = [y i1 , y i2 , ..., y in ] ∈ R n , where y ij is its estimation on the unit j's reference, j ∈ V . Under the LCA, we aim to design an algorithm satisfying where κ ij (0) is an initial value. We propose the Algorithm 3 as follows

Algorithm 3 With respecting power generation limits for a larger game
In the Algorithm 3, a ij is an entry of adjacent matrix A a (see Appendix A) and k > 0 is a constant. the difference between these estimate values of units i and k on the information of a unit j is defined as Additionally, q ij = a ij (−x j + c j ) and r ij = −a ijẋj for j ∈ P , and q ij = a ij c j and r ij = 0 for j ∈ L ∪ G ∪ T . We denote z i = [z i1 , ..., z in ] ⊺ ∈ R n as auxiliary estimation vector associated with player i and similarly define for Nj nẋ j for j ∈ P and d q ij (x) = 0 for j ∈ L ∪ G ∪ T . Thus, w ij z ij converges to w q ij z q ij = −(x j − c j ) for j ∈ P and to c j for j ∈ L ∪ G ∪ T , and d ij z ij converges to d q ij z q ij = −ẋ j for j ∈ P and to 0 otherwise, under initialization z ij (0) satisfying Lemma 6.2. By substituting these quasi-steady states into (21a), we have y r ij = (x j − c j ) + κ iP (0)e −γit for j ∈ P and y r ij = −c j + κ iP (0)e −γit for j ∈ L ∪ G ∪ T , where κ iP (0) = y r ij (0) − x r j (0) + c j and κ iP (0) = y r ij (0) − x r j (0) + c j . Then, satisfying the LCA requirement. We associate these equilibrium points (z * , µ * , w * , ν * , d * , η * , y * , λ * ) with the NE x * , and analyze the system (20)-(21) around these equilibriums. Let us obtain the following transformed system: We first derive a reduced equation forŷ ij asẏ r ij = −γ iŷ r ij + γ ixj +ẋ j if j ∈ P anḋ y r ij = γ iŷ r ij if j ∈ L ∪ G ∪ T . Thus,ŷ r ij =x j + θ r ij (0)e −γit for j ∈ P andŷ r ij = θ r ij (0)e −γit for j ∈ L ∪ G ∪ T , where θ r ij (0) =ŷ r ij (0) −x j (0). Similar to the analysis in Section 5, a reduced model is finally obtained aṡ with h i (t) = ρ ∑ j∈V,j≠i θ r ij (0)e −γit . We then construct a boundary-layer model from (23a) by applying the changes of variableŝ 1 .., (ŝ 2 n ) ⊺ ] ⊺ ∈ R n 2 and consider the time exchange expression τ = t to write the system (26) in a matrix form as In an analogous manner, one also can obtain a similar boundary-layer model for (23b)-(23c).

Proof. See Appendix B
In order to guarantee z ij ≠ 0 in (21a), we establish a condition as follows Lemma 6.2. If z ij is initialized at the neighborhood of their quasi-steady-state z q ij respectively, i.e., there exists a small constant > 0 such that max z ij (0) − z q ij ≤ , then z ij > 0 for ∀i, j ∈ V .

Proof. See Appendix B
Based on the analysis of boundary-layer model and reduced model, the final seeking result is stated as Theorem 6.3. Assume Slater condition, strict complementary slackness condition and φ i > ( P −3)ρ 2 are satisfied. There exists * > 0 such that 0 < i < * for ∀i ∈ V , if every unit in V adopts the proposed Algorithm 3 and their initializations satisfy Lemma 6.2 with λ li (0), λ ui (0) > 0, then each player's action eventually reaches the optimal power generation level locally in a distributed fashion.
Proof. See Appendix B

Remark 1. (Comparison between the Algorithm 2 and the Algorithm 3)
Although Algorithm 3 can be applied to solve more general games, it desires more computational resources than Algorithm 2, because each unit in Algorithm 3 maintains n-dimension real-time estimation vectors to exchange them with its neighbors, while in Algorithm 2 each unit needs only one variable, an estimation of average sum, to be exchanged. 7. Simulation in IEEE 118 bus system. IEEE 118 bus system provides a simple approximation of the American electric power system, comprising of 118 nodes, 91 loads, 54 generators, 9 transformers and 186 branches [9]. Regarding to this system, the size of each set is characterized by L = 54, T = 9, G = 17, P = 37. The local cost function for operating generator i is given as With Algorithm 1, estimation y i on the average of aggregate sum 1 n (x − ∑ i∈V c i ) is calculated in Fig. 1A, converging to the exact consensus value at −31.3568 (MW). In Fig. 1C Fig. 2A presents the estimation of 1 n (x − ∑ i∈V c i ) as −29.8339 (MW), which is achieved by all players in the game. Player units adopting this algorithm can seek their optimal power generation levels while respecting the given constraints as illustrated in Fig. 2B. For instance, in Fig. 2C Fig. 3A shows that all units' estimations on information of a player 4, which is x 4 − c 4 , exactly converge to the value value x * 4 − c 4 = −28.75 (MW). Moreover, Fig. 3B illustrates how the player 4 estimates information of other players in the game.
(b) Seeking NE by Algorithm 1.  8. Conclusion. In this paper, we developed a model for the SG based on the framework of IEEE 118 bus system. We showed that the problem of seeking the optimal set point for operating generators to achieve benefits selfishly turns out to be considered as an aggregative game in the interconnected network. To tackle the problem, we propose three distributed algorithms to reach that optimal set point, which is considered as a Nash Equilibrium of the game. The first one is designed when constraints are ignored, while the others take constraints into account. The simulation in the IEEE 118 bus system demonstrates the effectiveness of the proposed algorithms. There are some directions for future work. Firstly, other technical constraints associated to power flows and thermal phenomena may be considered. Secondly, we will extend our model for the network with uncertainties and switching topologies.
Appendix A. (Communication graph) The connected and undirected graph consists of n vertices, in which each vertex is presented for individual unit's smart meter,   .., G G , P 1 , ..., P P , T 1 , ..., T T } is the vertex set and E ⊆ V × V is the edge set. If (i, j) ∈ E then (j, i) ∈ E, and i, j are neighboring meters. We define N i = {j ∈ V (i, j) ∈ E} as a set of neighbors of the smart meter i. Associated with the graph, the adjacency matrix A a = [a ij ] ∈ R n×n is constructed according to interaction between meters, where a ij = 1 if (i, j) ∈ E, otherwise, a ij = 0. Moreover, the associated Laplacian matrix L = [l ij ] n×n is computed from L = A d − A a , where A d ∈ R n×n is a diagonal matrix with i-th element on the diagonal being ∑ j∈Ni a ij . Here, L1 n = 0 n and 1 ⊺ n L = 0 ⊺ n [20].

Appendix B.
Proof of Lemma 4.2: We firstly substitute quasi-steady statesẑ q i = − 1 n 1 ⊺ P x,ŵ q i = − 1 n 1 ⊺ P ẋ into slow dynamics (12), which yields a reduced solutionŷ r i = 1 n 1 ⊺ P x r + δ r i (0)e −γit and its derivative asẏ r i = 1 For convenience, we combine all exponentially decreasing functions to a collective vector r(t) = [e −γ1t , e −γ2t , ..., e −γ P t ] ⊺ ∈ R P , then writeŷ r i andẏ r i for all i ∈ P in a system aŝ  We then derive a dynamic equation forx i by substitutingŷ r i intoẋ i in (12) aṡ Finally, we collect all x r i for ∀i ∈ P to formẋ r and consider (28) to obtain the following reduced model. y r = −kBŷ r + kDdiag{δ r i (0)}r(t),ẋ r = −kBx r − kρndiag{δ r i (0)}r(t).

Proof of Theorem 4.3:
Here, the proof technique is based on singular perturbed algorithm in Section 2.1. We already mentioned that the conditions (B1) and (B2) are satisfied for the transformed system (11)- (12). Next, we are going to verify the condition (B3). Indeed, the boundary-layer model (13) is G.E.S, then for any bounded, andŷ i is bounded according to (12) as a result. It leads to that the right-hand side of system (11)-(12) and quasi-steady statesẑ q i ,μ q i ,ŵ q i ,ν q i are bounded. To prove the first partial derivatives of the right-hand side are bounded, we notice that the derivativesż i ,μ i ,ẇ i ,ν i ,ẏ i are bounded, thus partial derivatives of the right-hand side according to variable t is bounded, so are partial derivatives according to other variables due to a chain rule, i.e., ∂ẑi ∂xi = ∂ẑi ∂t ∂t ∂xi , which is bounded.
Also, the first partial derivatives of quasi-steady states are bounded. Similarly, one can conclude that second partial derivatives of the right-hand side and quasi-steady states are bounded, which satisfy condition (B3). To conclude the proof, we utilize the result in the globally exponential stable property of boundary-layer model and ∂xi∂xj ∈ R P × P as strictly diagonally dominant matrices with positive diagonal entries as well as symmetric matrices, and thus positive definite. Therefore, H(x, α) and f (x i ,x, α) are strictly convex in x, and the optimization problem (5) and (14) are equivalent. ◻ Proof of Lemma 5.2: Here, we begin with the unforced system: whereσ r is bounded andσ r globally asymptotically converges to the origin ([18]-Lemma 1). In addition, the systemẏ r = f y (ŷ r ,σ r ) = −kBŷ r + k n 1 P × P (λ r l −λ r u ), withλ r l −λ r u considered as an input and −kB is Hurwitz matrix, is input-to-state stable system since there exists a class KL function β and a class K function υ so that its solution satisfies: ŷ r (t) ≤ β( ŷ r (t 0 ) , t − t 0 ) + υ sup t0≤η≤t λ r l (η) −λ r u (η) for any initial stateŷ(t 0 ) and any bounded inputλ r l −λ r u . Then, since lim t→+∞ (λ r l − λ r u ) = 0, the systemẏ r = f y (ŷ r ,σ r ) is G.A.S [11], so is (30), meaning that the trajectories of the unforced system (30) is proved to globally converge to the origin.
Because it can be directly verified from (18) thatf y (ŷ r ,σ r , t) andf σ (σ r , t) are bounded for boundedŷ r andσ r as r(t) and Ξ(t) are exponentially decaying terms, in order to apply Limiting equation theorem in Section 2.2 to prove the reduced model (18) be G.A.S., it remains to show the trajectoriesŷ r andσ r of the reduced model (18) are bounded.