A radial distribution network is represented by a directed graph G p . Except for the substation node having a constant voltage, buses in the distribution networks are denoted by the set N , among which the buses with storage are collected in the set S , while all the other buses are collected in the set ℒ = N \ S . Note that the buses in ℒ involve loads and renewables. The set of line segments connecting buses in N are denoted by the edge set Ε p = { ( i , j ) ⊂ N × N } . For every line segment ( i , j ) ∈ Ε p , the DistFlow model is established as follows: P i j − ∑ k ∈ N j P j k = − P j n , (1a) Q i j − ∑ k ∈ N j Q j k = − Q j n , (1b) V i − V j = r i j l P i j + x i j l Q i j , (1c) where P i j , Q i j denote the active and reactive power flows on line ( i , j ) ; the net injections of the bus j are P j n = P j g − P j l ,     Q j n = Q j g − Q j l , where P j g , Q j g denote active/reactive powers from DERs and P j l , Q j l denote the powers of load consumptions. The power flows from the substation to the end buses lead to voltage drops as V i − V j described in (1c), where r i j l , x i j l are resistance and reactance of the distribution lines, N j denotes the set of neighbors of note j . The DistFlow model in (1) can be represented in a more compact form as (Zhu and Liu, 2016) follows: V = R l P n + X l Q n + V 0 , (2a) R l = M − T D r M − 1 ,   X l = M − T D x M − 1 . (2b)

In Equation 2, V is the vector of bus voltages and denoted as V = ( V j ) j ∈ N ; similarly, P n = ( P j ) j ∈ N ,     Q n = ( Q j ) j ∈ N are vectors of net bus injections, and V 0 = v 0 1 | N | where v 0 is the voltage magnitude of node 0 and 1 | N | denotes a vector of all one with dimension | N | , where | N | is the cardinality of N . The matrix M is the incidence matrix of the graph ( N , ℰ p ) , and D r = diag ( r i j l ) i j ∈ ℰ p , D x = diag ( x i j l ) i j ∈ ℰ p are diagonal matrices of the line parameters.

The cyber network communicating distributed storage devices are represented as a graph G c = { S , ℰ c } in this study, where the graph edges in ℰ c represent the communication lines among storage agents in S . The connectivity of the graph is described by the adjacency matrix A c = ( a i j ) i , j ∈ S where a i j = 1 for all j ∈ S i and 0 otherwise, where S i is the neighboring set of node i . We assume that the cyber graph is undirected and connected.

The optimization problem is formulated for the optimal storage coordination of a VSP located in a radial distribution network. The VSP is controlled to deliver a given amount of active power while regulating the local network voltage using its reactive power control capability. Considering every time slot t in a time horizon T and ∀ i ∈ S , we have min . ( P i , t b , Q i , t b ) ∑ i = 1 | S | { α i , t P ( P i b ) 2 + γ i , t P | P i b | + α i , t Q ( Q i b ) 2 + γ i , t Q | Q i b | } , (3a) V min 1 | S | ≤ ∑ i = 1 | S | { R i l ( P i , t b + P i , t l ) + X i l ( Q i , t b + Q i , t l ) } + ∑ j = 1 | ℒ | { R j l P j , t n + X j l Q j , t n } + V 0 ≤ V max 1 | S | , (3b) ∑ i = 1 | S | P i , t b = P t A , (3c) ( P i , t b ) 2 + ( Q i , t b ) 2 ≤ ( S i b ) 2 , (3d) P i b , m i n ≤ P i , t b ≤ P i b , m a x , (3e) S o C m i n ≤ S o C i , t b = S o C i , t − 1 b + η P i , t b C A i b ≤ S o C m a x ,       ∀ i ∈ S . (3f)

In the optimization problem (3), P i , t b , Q i , t b are active and reactive power outputs of the i th storage unit at the current time step. We assign a quadratic cost for every storage unit that is described by the parameters α i , t P ,   γ i , t P ,   α i , t Q ,   γ i , t Q as in (Eq. 3a), and the absolute values indicate that both charging and discharging yield costs are related to storage usage. In (Eq. 3b), R i l = ( R i j l ) j ∈ S ,   X i l = ( X i j l ) j ∈ S ∈ R | S | are column vectors of R l ,   X l corresponding to the i th storage response on other storage buses, which can be understood as the bus voltage sensitivities with respect to the power injections at the i th buses. Here, we separately represent the power injections on the storage buses and load buses in ℒ , whose net powers P j , t n ,   Q j , t n combine the effects of renewables and loads as in (Eq. 3b); V m i n ,   V m a x are the bus voltage limits. The VSP is optimized to provide the aggregated power P t A as in (Eq. 3c). We mainly focus on inverter-based storage technologies, and the local constraints (3d)—(3e) characterize the limits of storage power capacity where S i b is the apparent power limit of the power converter. Eq. 3f imposes the limits on storage energy capacity based on a general storage model, in which η is the coefficient for charging/discharging, and C A i b is the storage energy capacity. Note that (Eq. 3f) can be incorporated into (Eq. 3e) with time-varying power limits depending on the storage SoC (Dall’Anese et al., 2018); for example, storage with a full charge has P i , t b , m a x = 0 and P i , t b , m i n = − S i b . To facilitate the algorithm implementation, the quadratic constraint (3d) is further linearized as follows: − S i b ≤ cos ( τ π κ ) P i , t b + sin ( τ π κ ) Q i , t b ≤ S i b ,       τ = 1 , … , κ , (4) where the accuracy loss is 1.5% when κ = 8 (Jabr, 2017). Moreover, considering that only the power injections on storage buses are known by agents, the network voltage constraint (3b) is rewritten in an increment form of storage power, as follows: V m i n 1 | S | ≤ ∑ i = 1 | S | { R i l ( P i , t − 1 b + Δ P i , t b + P i , t l ) + X i l ( Q i , t − 1 b + Δ Q i , t b + Q i , t l ) } + ∑ j = 1 | ℒ | { R j l P j , t n + X j l Q j , t n } + V 0 ≤ V m a x 1 | S | . (5)

It is assumed that the algorithm is fast enough such that the powers of loads and renewables are nearly constant within each time slot. From (Eq. 5) we denote V ˜ t = ∑ i = 1 | S | { R i l ( P i , t − 1 b + P i , t l ) + X i l ( Q i , t − 1 b + Q i , t l ) } + ∑ j = 1 | ℒ | { R j l P j , t n + X j l Q j , t n } + V 0 , (6) which is the vector of voltage magnitudes of storage buses before the agent actions and can be measured at the beginning of the current time step t . We denote x i , t = ( P i , t b , Q i , t b ) T = x i , t − 1 + x i , t , and then the optimization problem (3) is arranged in the canonical form with the incremental variables x i , t , considering storage outputs at t − 1 are known by each agent. min Δ x i , t , Δ r i , t ≥ 0 ∀ i ∈ S ∑ i = 1 | S | { Δ x i , t T Ω i , 0 Δ x i , t + Ω i , t − 1 Δ x i , t + ‖ Γ i Δ x i , t + Γ i x i , t − 1 ‖ } , (7a) ∑ i = 1 | S | ( [ A i E ] Δ x i , t + [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] ) ≤ 0 , (7b) C Δ x i , t + Δ r i , t − d i = 0 ,       ∀ i ∈ S , (7c) where Ω i , 0 = d i a g ( α i , t P , α i , t Q ) ,   Ω i , t − 1 = 2 x i , t − 1 T Ω i , 0 ,   Γ i = d i a g ( γ i , t P , γ i , t Q ) , A i = [ R i l , X i l − R i l , − X i l ] ,   b i , t − 1 = [ V ˜ i , t 1 i − V ˜ i , t 1 i ] ,   b 0 = [ V max 1 | S | − V min 1 | S | ] , E = [ 1 0 − 1 0 ] ,   e i , t − 1 = [ E x i , t − 1 − E x i , t − 1 ] ,   e 0 = [ P t A − P t A ] , C s = [ cos ( 1 π κ ) sin ( 1 π κ ) cos ( 2 π κ ) sin ( 2 π κ ) ⋮ ⋮ cos ( κ π κ ) sin ( κ π κ ) ] , C p = [ 1,0 ] , C = [ C s − C s C p − C p ] , d i = [ S i b 1 κ − C s x i , t − 1 S i b 1 κ + C s x i , t − 1 P i , t max − C p x i , t − 1 − P i , t min + C p x i , t − 1 ] .

The vector 1 i has the i th entity equal to 1 and 0 for the others. The 1-norm ‖ ■ ‖ 1 in the objective function (7a) accounts for storage costs in both charging and discharging modes; the coupled equality constraints 3(c) are represented in inequality, which with (Eq. 3b) is incorporated into (7b); (7c) are local constraints for storage capacity limits while Δ r i , t ≥ 0 is a slake variable.

Consider the following optimization problem: min . x , r h ( x ) + g ( r ) , (8a) A x + B r = c . (8b)

For the augmented Lagrangian of (8), as L ρ ( x , r , y ) = h ( x ) + g ( r ) + y T ( A x + B r − c ) + ( ρ / 2 ) || A x + B r − c 2 2 , ADMM consists of the following iterations: x k + 1 = arg min . x L ρ ( x , r k , y k ) , (9a) r k + 1 = arg min . r L ρ ( x k + 1 , r , y k ) , (9b) y k + 1 = y k + ρ ( A x k + 1 + B r k + 1 − c ) . (9c)

The updates of the primal variables are carried out jointly with respect to the augmented Lagrangian in a sequential fashion, while the dual update takes on a gradient step with coefficient ρ . ADMM brings robustness and yields convergence without the restrictive assumptions for primal-dual and dual ascent algorithms (Boyd et al., 2010).

In the classical ADMM, a central coordinator is usually needed to collect the primal variables in (Eq. 9c), and then the updated dual variable is broadcasted to every agent. To make the ADMM completely distributed, a consensus ADMM is proposed by Mateos et al. (2010), which is further developed to solve the problem with coupled constraints via its Lagrange dual problem (Chang et al., 2015;Chang, 2016). For an optimization problem in the form of (7), in which we make the generalization that Δ x i , t ∈ R n ,   Δ r i , t ∈ R m and the other matrices have compatible sizes. The objection function of (7) can be separated for individual agents and denoted as f i , t ( Δ x i , t ) = g i , t ( Δ x i , t ) + h i , t ( Δ x i , t ) , where g i , t ( Δ x i , t ) and h i , t ( Δ x i , t ) are the smooth and non-smooth parts, respectively. The Lagrange dual problem of (7) is considered and arranged in a consensus form, as follows: min . y i , t ≥ 0 , z i , t   max . Δ x i , t , Δ r i ≥ 0 i ∈ S ∑ i = 1 S { − f i ( Δ x i , t ) − y i , t T ( [ A i E ] Δ x i , t + [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] ) − z i , t T ( C Δ x i , t + Δ r i , t − d i ) } , (10a) y i , t = t i j , y j , t = t i j , ∀ j ∈ S i , (10b) z i , t = l i , i ε S . (10c) where y i , t ,   z i , t are dual variables associated with the coupled constraints and local constraints. Since the dual variables of the coupled constraints are kept by individual agents, constraint (10b) imposes their consensus, in which t i j is an auxiliary variable. Constraint (10c) is a dummy constraint to make a strongly convex subproblem facilitating the following derivation (Chang, 2016). The augmented Lagrangian of (10) is as follows: L c ( y i , z i , λ + , i j , λ − , i j , ν i ) = ∑ i = 1 S { φ i ( y i , z i ) + 1 | S | y i T ( [ b i , t − 1 e i , t − 1 ]   − [ b 0 / | S | e 0 / | S | ] ) + z i T d i   + ∑ j ∈ S i [ λ + , i j T ( y i − t i j ) + λ − , i j T ( y j − t i j ) ]   + ν i T ( z i − l i ) + τ 2 ‖ z i − l i ‖ 2 2   + σ 2 ∑ j ∈ S i ( ‖ y i − t i j ‖ 2 2 + ‖ y j − t i j ‖ 2 2 ) } . (11) where we drop the time stamp t to simplify notions, and φ i ( y i , z i ) = max . Δ x i , Δ r i ≥ 0 { − f i ( Δ x i ) − y i T [ A i E ] x i , t − z i T ( C i Δ x i + Δ r i ) } , and λ + , i j T ,   λ − , i j T , ν i are dual variables associated with the consensus and dummy constraints, τ ,   σ are positive coefficients of the Lagrange multipliers. Applying the ADMM procedure for (11), we have ( y i , z i ) = arg   min . y i ≥ 0 , z i { φ i ( y i , z i ) + 1 | S | y i T ( [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] )   + z i T d i + ∑ j ∈ S i [ λ + , i j T ( y i − t i j ) + λ − , i j T ( y j − t i j ) ]   + σ 2 ∑ j ∈ S i ( ‖ y i − t i j ‖ 2 2 + ‖ y j − t i j ‖ 2 2 )   + τ 2 ‖ z i − l i k − 1 + ν i k − 1 τ ‖ 2 2 } , (12a) t i j k = arg   min . t i j { ∑ j ∈ S i ( ‖ y i k − t i j + λ + , i j k − 1 σ 1 ‖ 2 2 + ‖ y j k − t i j + λ − , i j k − 1 σ 1 ‖ 2 2 ) } , (12b) l i k = arg min . l i { ‖ z i k − l i + ν i k − 1 τ ‖ 2 2 } , (12c) ν i k = ν i k − 1 + τ ( z i k − l i k ) , (12d) λ + , i j k = λ + , i j k − 1 + σ 1 ( y i k − t i j k ) ,       λ − , j i k = λ − , j i k − 1 + σ 1 ( y i k − t j i k ) . (12e)

Deriving solutions of (12a)–(12c) and combining with (12d)–(12e) gives t i j k = t j i k = ( y i k + y j k ) / 2 , z i k = l i k ,   ν i k = 0 . So, (12) can be simplified as follows: ( y i , z i ) = arg min . y i ≥ 0 , z i { φ i ( y i , z i ) + 1 | S | y i T ( [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] )   + z i T d i + y i T ∑ j ∈ S i ( λ + , i j k − 1 + λ − , i j k − 1 )   + ∑ j ∈ S i ( σ ‖ y i − y i k − 1 + y j k − 1 2 ‖ 2 2 ) + τ 2 ‖ z i − z i k − 1 ‖ 2 2 } , (13a) λ + , i j k = λ + , i j k − 1 + σ ( y i k − y i k + y j k 2 ) ,   λ − , j i k = λ − , j i k − 1 + σ ( y i k − y i k + y j k 2 ) . (13b)

Using the minmax theorem and the strong convexity of (13a) (Chang, 2016), y i , z i can be conveniently obtained by expanding φ i and combining the linear and quadratic terms, which gives ( Δ x i , t k , Δ r i , t k ) = arg min . Δ x i , t , Δ r i , t ≥ 0 { f i , t ( Δ x i , t ) + σ 4 | S i | ‖ 1 σ ( [ A i E ] Δ x i , t + [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] ) − p i , t k − 1 σ + ∑ j ∈ S i ( y i , t k − 1 + y j , t k − 1 ) ‖ 2 2 + 1 2 τ ‖ C Δ x i , t + Δ r i , t − d i + τ z i , t k − 1 ‖ 2 2 } , (14a) y i , t k = 1 2 | S i | [ ∑ j ∈ S i ( y i k − 1 + y j k − 1 ) − 1 σ p i , t k − 1 + 1 σ ( [ A i E ] Δ x i , t + [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] ) ] + , (14b) z i , t k = z i , t k − 1 + 1 τ ( C Δ x i , t k + Δ r i , t k − d i ) , (14c) p i , t k = p i , t k − 1 + c ∑ j ∈ S i ( y i , t k + y j , t k ) , (14d) where p i k = ∑ j ∈ S i ( λ + , i j k + λ − , j i k ) and [ ■ ] + = max ( 0 , ■ ) . The iterative minimization of primal variables in (14a) does not need to be very accurate since it is an intermediate step in ADMM procedures (Mateos et al., 2010). In this regard, the inexact ADMM approximates (14a) by its first-order Taylor expansion at ( Δ x i , t k − 1 , Δ r i , t k − 1 ) (Chang, 2016). Denote the smooth part of (14a) as follows: g ˜ i , t ( Δ x i , t ) = Δ x i , t T Ω i , 0 Δ x i , t + Ω i , t − 1 Δ x i , t           + σ 4 | S i | ‖ 1 σ ( [ A i E ] Δ x i , t + [ b i , t − 1 e i , t − 1 ] − [ b 0 / | S | e 0 / | S | ] ) − 1 σ p i , t k − 1 + ∑ j ∈ S i ( y i , t k − 1 + y j , t k − 1 ) ‖ 2 2 + 1 2 τ ‖ C Δ x i , t + Δ r i , t − d i + τ z i , t k − 1 ‖ 2 2 . (15)

To derive an analytical solution, (14a) is simplified as follows: ( Δ x i , t k , Δ r i , t k ) = arg min . Δ x i , t , Δ r i , t ≥ 0 { h i , t ( Δ x i , t ) + ( ∇ x g ˜ i , t k − 1 ) T ( Δ x i , t   − Δ x i , t k − 1 ) + ( ∇ r g ˜ i , t k − 1 ) T ( Δ r i , t − Δ r i , t k − 1 )   + β i 2 ‖ Δ x i , t − Δ x i , t k − 1 ‖ 2 2 + β i 2 ‖ Δ r i , t − Δ r i , t k − 1 ‖ 2 2 } , (16) where β i > 0 is some penalty coefficient. Expanding the gradients ∇ x g ˜ i , t k − 1 and ∇ r g ˜ i , t k − 1 gives Δ x i , t k = arg min . Δ x i , t { ‖ Γ i Δ x i , t + Γ i x i , t − 1 ‖ 1 + β i 2 ‖ Δ x i , t − ( Δ x i , t k − 1 − 1 β i ∇ x g ˜ i , t k − 1 ) ‖ 2 2 } , (17a) Δ r i , t k = [ ( 1 − 1 β i τ ) Δ r i , t k − 1 + 1 β i τ ( d i − C Δ x i , t k − 1 − τ z i , t k − 1 ) ] + . (17b)

To solve (17a) with non-smooth 1-norm terms, we employ the proximal operator for each coordinate of Δ x i , t . Denoting the j th component of Δ x i , t k as Δ x i , t k ( j ) , we have Δ x i , t k ( j ) = arg min . Δ x i , t ( j ) { | Γ i ( j ) Δ x i , t ( j ) + Γ i ( j ) x i , t − 1 ( j ) | + β i 2 [ Δ x i , t ( j ) − ( Δ x i , t k − 1 ( j ) − 1 β i ∇ x g ˜ i , t k − 1 ( j ) ) ] 2 } , (18) where Γ i ( j ) is the j th diagonal element of Γ i and ∇ x g ˜ i , t k − 1 ( j ) is the j th component of the gradient vector. Since h i , t ( Δ x i , t ) = ∑ j = 1 n h i , t j ( Δ x i , t ( j ) ) = ∑ j = 1 n | Γ i ( j ) Δ x i , t ( j ) + Γ i ( j ) x i , t − 1 ( j ) | is separatable, and notice that (Eq. 18) is a proximal operator for a scaler function h i , t j ( Δ x i , t ( j ) ) , which is precomposed by an affined form. Thus, we rearrange (18) as follows: Δ x i , t k ( j ) = arg min . Δ x i , t ( j ) { h i , t j ( Δ x i , t ( j ) )   + [ Δ x i , t ( j ) − ( Δ x i , t k − 1 ( j ) − 1 β i ∇ x g ˜ i , t k − 1 ( j ) ) ] 2 } = p r o x . h i , t j / β i ( Δ x i , t k − 1 ( j ) − 1 β i ∇ x g ˜ i , t k − 1 ( j ) ) = 1 Γ i ( j ) { p r o x Γ i 2 ( j ) | ■ | / β i [ Γ i ( j ) ( Δ x i , t k − 1 ( j )   − 1 β i ∇ x g ˜ i , t k − 1 ( j ) ) + Γ i ( j ) x i , t − 1 ( j ) ] − Γ i ( j ) x i , t − 1 ( j ) } , (19) where the last equation uses the scaling and translation property of the proximal operator (denoted as prox.) for scalar functions (Beck, 2017). The proximal operator of the absolute value is called soft thresholding (Boyd et al., 2010), and thus, (19) can be computed analytically as follows: Δ x i , t k ( j ) = 1 Γ i ( j ) { [ Γ i ( j ) ( Δ x i , t k − 1 ( j ) − ∇ x g ˜ i , t k − 1 ( j ) β i ) + Γ i ( j ) x i , t − 1 ( j ) − Γ i 2 ( j ) β i ] + − [ − Γ i ( j ) ( Δ x i , t k − 1 ( j ) − 1 β i ∇ x g ˜ i , t k − 1 ( j ) ) − Γ i ( j ) x i , t − 1 ( j ) − Γ i 2 ( j ) β i ] + − Γ i ( j ) x i , t − 1 ( j ) } . (20)

From (Eq. 20), we recover the optimal VSP control as defined in the System Modeling and Problem Formulation section. This gives the control commands for storage agents: Δ P i , t k = 1 γ i , t P { [ γ i , t P ( Δ P i , t k − 1 − 1 β i ∇ P g ˜ i , t k − 1 ) + γ i , t P P i , t − 1 − γ i , t P   2 β i ] + − [ − γ i , t P ( Δ P i , t k − 1 − 1 β i ∇ P g ˜ i , t k − 1 ) − γ i , t P P i , t − 1 − γ i , t P   2 β i ] + − γ i , t P P i , t − 1 } , (21a) Δ Q i , t k = 1 γ i , t Q { [ γ i , t Q ( Δ Q i , t k − 1 − 1 β i ∇ Q g ˜ i , t k − 1 ) + γ i , t Q Q i , t − 1 − γ i , t Q   2 β i ] + − [ − γ i , t Q ( Δ Q i , t k − 1 − 1 β i ∇ Q g ˜ i , t k − 1 ) − γ i , t Q Q i , t − 1 − γ i , t Q   2 β i ] + − γ i , t Q Q i , t − 1 } . (21b)

Combining (21) with the primal and dual updates in (17b), (14b)–(14d) makes up the control law of storage devices in a VSP. Assuming a zero-duality gap, and the primal and dual optimal is attainable, the only requirement for the algorithm convergence is that the penalty coefficient β i should be larger than some constant determined by the modulus and Lipschitz constant of the smooth part of the objection function g i , t (Chang, 2016).

The simulations use real PV generation data sampled at 5 min, which are collected from a 75 MW solar power plant in Colorado, United States (Solar Power Data for Integration). The PV generation data from different days are used (Solar Power Data for Integration Studies, 2006) for 25 PV systems covering 24 h, and their profiles are shown in Figure 2. The PV powers are scaled-down by 150 times to suit the capacities of the distribution networks. The buses where the storage and PVs are located are randomly selected, while the loads are assumed to be constant in the simulations. The VSP dispatch commands are each randomly generated in 10 min. The voltage limits for cases 1 and 2 are 0.95–1.05 p.u., and 0.99–1.01 p.u. for case 3. Figure 3 shows the voltage profile of the IEEE 33 and IEEE 69 bus systems with integrated PV systems. This indicates that without storage response, a significant voltage increase is observed during peak PV generation hours.

Table 1 gives the storage and control parameters in the simulations. The storage agents are labeled by the same bus number where they are located. We assume a ring topology for storage communication in the VSP cyber network, that id each storage agent communicates to its 2 most adjacent neighbors in the undirect graph. The sampling time of the algorithm is set to 1 min, which gives enough time for algorithm convergence at each time step. The storage power capacities and costs are randomly selected from the ranges provided in Table 1. We assume that all the storages have enough energy in the considered duration. The condition for algorithm termination is the accuracy of the objection function ( obj k − obj ∗ ) / obj ∗ ≤ 1 e − 5 , where obj ∗ is the cost of a centralized optimization for the same problem.

Case I: The IEEE 33-bus system is rated at 12.7 kV with 3.75 MW/2.3MVar loads in total. The PV generations are selected from the data profile and integrated into the network at random buses. Figures 4A, B compares the VSP aggregated powers and costs under the proposed distributed control and a centralized optimization using YALMIP. This illustrates that the proposed VSP control can track the given power reference in a timely fashion and yield almost identical storage costs with respect to a centralized optimization result. The voltage profiles of the system are given in Figures 4C,D. Again, the proposed distributed control, i.e., Figure 4D generates similar results as the centralized scenario in Figure 4C, both of which successfully restrict the voltage within 0.95–1.01p.u., compared to the significant voltage variation observed in Figure 2A. Generally, the computational time of the distributed approach scales linearly with respect to the dimension of the decision variables, while the centralized approach has a high-order polynomial relation between the computational time and the dimension.

Case II: The IEEE 69 bus system is used in this study case. The radial distribution network is rated at 12.7 kV. The total loads are 3.8 MW and 2.7 MVar. Again, the simulation results in Figure 5 indicate that the proposed distributed control accomplishes nearly identical results as a centralized optimization in VSP dispatch and voltage regulations, while the distributed control has fewer requirements on the communication capability of the cyber network.

Case III: The data of the 119-bus system is obtained from Zhang et al. (2007) and the diagram is shown in Figure 6. The nominal voltage of the system is 11 kV, and all the initial tie switches are open to making a radial distribution network. The total power loads are 22.7 MW and 17.0 MVar. In Figure 7, we present the active powers and reactive powers of the total 40 storage units in the VSP. The first column gives the results from a centralized optimization using CPLEX, while the second column is the result using the proposed PDC-ADMM method. It can be observed that the centralized control maintains the storage power references in each dispatch duration, while the proposed distributed control can gradually approach the optimal dispatch using feedback from the previous time step. Figure 8 validates the performance of the proposed VSP control in terms of the power reference tracking and voltage regulations similar to the results in IEEE 33-bus and IEEE 69-bus systems. Note that the proposed control framework generalizes the functionalities of storage solely by providing voltage supports in distribution networks.

The centralized optimization can provide the global optimal regarding the storage dispatch in the convexified optimization problem. It provides a lower bound on the cost of storage coordination. The proposed distributed control framework uses an inexact approach to reduce the computational efforts of each storage agent. So, the accuracy of the optimization is compromised in exchange for the speed of storage dispatch. In addition, for the sake of practical implementation, we set a fixed value for the maximum number of iterations. So, the distributed control approach might output the calculated storage setpoints even when a suboptimal solution is obtained. These are the reasons that cause slight differences in the cost and storage dispatch between the centralized and distributed approach, but the control performance regarding the reference tracking and voltage regulation is almost in the two control scenarios.