Distributed Optimization-based Electric Vehicle Charging and Discharging in Unbalanced Distribution Grids

—The worldwide proliferation of Electric Vehicles (EVs) is accelerating the need for distributed optimization approaches that are both scalable and computationally efﬁcient, to enable coordinated EV charging and discharging. This paper proposes a distributed optimization-based algorithm called Dis-Net-EVCD for network-aware EV charging and discharging in unbalanced distribution grids, incorporating grid voltage con- straints and EV customer economics. The Alternating Direction Method of Multipliers (ADMM) underpins the development of our Dis-Net-EVCD algorithm, wherein EVs determine their charge-discharge proﬁles locally via peer-to-peer communication. A receding horizon implementation of Dis-Net-EVCD supports near-real-time EV charging and discharging, while accommodat- ing unexpected EV arrivals and departures. Numerical simulations carried out on the IEEE 13 node test feeder demonstrate that EV customers implementing Dis-Net-EVCD yield a total operational cost reduction of 78% compared to uncoordinated EV charging, while improving grid voltage regulation services and satisfying EV charge requirements ahead of departure. Moreover, in simulation, Dis-Net-EVCD is shown to be approximately 60 times computationally faster than its centralized counterpart.

The studies in [7,8] propose distributed EV charging schemes based on the shrunken primal-dual subgradient algorithm to fill the load-curve valley while operating within compliant steady-state voltage constraints. Alternating Direction Method of Multipliers (ADMM)-based distributed EV charging schemes are proposed in [9][10][11] to avoid system voltages dropping below the permissible limits. Other distributed EV charging schemes focused on reducing voltage fluctuations in a distribution grid include droop control [12], game theory [13] and multi-agent system modelling [14]. While the aforementioned literature is limited to EV charge-only operations, the distributed approaches in [15,16] incorporate Vehicle-to-Grid (V2G) operation to further improve voltage regulation.
In practice, three-phase distribution grids are typically unbalanced in both load and impedance across phases, especially when designed with single phase and two phase laterals. With the exception of [10], the distributed EV charging schemes in [7][8][9][10][11][12][13][14][15][16] incorporate single-phase representations of threephase distribution grids. However, single-phase representations of the power flow equations cannot capture the steady-state voltage variability across phases in an unbalanced distribution grid. Interestingly, the distributed EV charging scheme in [10] minimizes power supply costs subject to regulating unbalanced steady-state voltages. Similar to [10], we seek to develop a distributed approach that yields minimal costs for EV charging while regulating quasi-steady-state voltages across phases. Beyond [10], we seek to include V2G operations and remove the requirement for a central coordinator -reducing the risk of communication-based single point of failures.
In this paper, we propose an algorithm for Distributed and Network-aware EV Charging and Discharging, referred to as Dis-Net-EVCD. Specifically, Dis-Net-EVCD minimizes EV customer operational costs while maintaining the voltages in an unbalanced distribution grid within prescribed quasi-steadystate limits. We consider EV customer operational costs associated with: (1) purchasing (or otherwise being compensated for delivering) electricity on a Time-of-Use (ToU) net-metering tariff; and (2) battery degradation due to frequent charging and discharging. Dis-Net-EVCD also accommodates individual EV preferences and battery constraints, including EV charge-discharge rate limits, EV battery state-of-health thresholds, and customer-specified charging requirements. Central to Dis-Net-EVCD is the Dual Consensus -ADMM (DC-ADMM [17]) decomposition method, in which a centralized optimization problem is decomposed into a set of tractable subproblems that are solved by EVs in parallel by exchanging information with neighboring EVs.
The contributions of this paper are summarized as follows: 1) By leveraging from DC-ADMM -combining ADMM, duality theory and consensus theory -we propose a distributed and network-aware algorithm for EV charging and discharging called Dis-Net-EVCD. The objective is to minimize operational costs accrued to EV customers while regulating voltage magnitudes in unbalanced distribution grids to within permissible limits. That is, Dis-Net-EVCD carefully balances reductions in EV operational costs against improvements in supply voltages across the distribution grid, and thereby benefits both EV customers and electrical distributors. Importantly, Dis-Net-EVCD supports peer-to-peer communication between EVs, eliminating the need for a central coordinator. 2) We prove that Dis-Net-EVCD asymptotically converges to the globally optimal solution, where EVs seek to reduce their operational costs and fulfill their heterogeneous charge requirements -all while complying with EV battery constraints and grid voltage constraints.

3) We propose a receding-horizon implementation of
Dis-Net-EVCD that manages non-deterministic EV arrivals and departures. The remainder of the paper is organized as follows. Section II introduces a residential EV energy system. Section III formulates a centralized optimization problem for networkaware EV charging and discharging, followed by the proposed distributed algorithm in Section IV. Section V presents our numerical simulations and Section VI concludes the paper.

Notation
Let R denote -dimensional vectors of real numbers and let R ≥0 denote -dimensional vectors with all non-negative elements, where, as usual, R 1 = R. Let A denote the transpose of a vector or matrix A. Let 0 ∈ R denote the all-zero column vector of length and let 1 ∈ R denote the all-ones column vector of length . Let I ∈ R × denote the -by-identity matrix. For two matrices A = [a ij ] ∈ R m×n and B = [b ij ] ∈ R p×q , let A ⊕ B := diag(A, B) = [c ij ] ∈ R (m+p)×(n+q) denote the matrix direct sum, satisfying c ij = a ij when i ≤ m and j ≤ n, c ij = b (i−m)(j−n) when i > m and j > n, and c ij = 0 elsewhere. Let

II. PRELIMINARIES
A. Residential EV Energy System (REES) Fig. 1 illustrates the topology of a Residential EV Energy System (REES) of a single customer, consisting of an EV and residential load situated behind the Point of Common Coupling (PCC). The residential load represents either the energy consumed by the customer, or any excess generation produced by co-located rooftop solar. The set of REESs are indexed by N := {1, . . . , n, . . . , N }, where customer n ∈ N charges (and discharges) EV n .
Consider a planning time-horizon [0, T ∆] consisting of T time-intervals of length ∆ (in hours). Let S := {1, . . . , t, . . . , T }, where a time-interval is denoted by (t − 1)∆, t∆ . In this paper, we consider a 24-hour time-horizon Fig. 1. Residential EV energy system (REES) illustrating the electrical path from customer n ∈ N to the wider distribution grid. Bi-directional smart meter M measures and records power flow gn(t) for the purpose of electricity billing and compensation. For each n ∈ N, the power balance equation gn(t) = ln(t) + xn(t) holds for all t ∈ S, where arrows associated with xn(t), ln(t) and gn(t) represent the assumed direction of positive power flow.
with ∆ = 0.5 and T = 48. Other choices are certainly possible, subject to commensurability of T and ∆.
The power flows in Fig. 1 are as follows. The measured power (in kW) through meter M is denoted by g n (t) ∈ R ≥0 (or g n (t) ∈ R <0 ), which is the power delivered from (to) the grid to (from) customer n ∈ N over the time-interval (t − 1)∆, t∆ . Accordingly, the grid power profile of customer n ∈ N over [0, T ∆] is defined as g n := [g n (1), . . . , g n (t), . . . , g n (T )] ∈ R T . Throughout, EVindices are sub-scripted and time-indices are parenthesized.
The non-EV residential power consumption (or excess renewable power generation) (in kW) of customer n ∈ N over the time-interval (t−1)∆, t∆ is represented by l n (t) ∈ R ≥0 (or l n (t) ∈ R <0 ) for all t ∈ S. In more detail, l n (t) < 0 when the rooftop solar generation is in excess of the non-EV power consumption, and l n (t)> 0 when power is delivered from the grid (or EV) to the residential load. Accordingly, the dayahead non-EV power profile of customer n ∈ N over [0, T ∆] is defined by l n := [l n (1), . . . , l n (t), . . . , l n (T )] ∈ R T . Likewise, the power (in kW) delivered to (from) EV n over the time-interval (t − 1)∆, t∆ is denoted by x n (t) ∈ R ≥0 (or x n (t) ∈ R <0 ) for all t ∈ S. Accordingly, the chargedischarge profile of EV n over [0, T ∆] is defined by x n := [x n (1), . . . , x n (t), . . . , x n (T )] ∈ R T , where x n (t) > 0 when the EV charges and x n (t) < 0 when the EV discharges.
The EV of each customer n ∈ N, EV n , is associated with a set of design parameters that include the arrival time index a n ∈ S, intended departure time index d n ∈ S, battery capacity c n (in kWh), initial State-of-Charge (SoC)σ n (in kWh), target SoC σ * n (in kWh), minimum SoC σ n (in kWh), maximum SoC σ n (in kWh), maximum charge rate x n (in kW), maximum discharge rate x n (in kW), charge efficiency µ n (0 ≤ µ n ≤ 1), and discharge efficiency µ n (µ n ≥ 1). Clearly, if EV n is charge-only, then x n = 0. Upon arrival at home, the EV owner must specify an expected departure time index, d n , and a target SoC that is feasible within the expected departure time, i.e., the condition σ * n ≤σ n + (d n − a n )x n µ n ∆ must hold [2,18]. The collective charging specifications for EV n are represented by Λ n := {a n , d n , c n ,σ n , σ * n , σ n , σ n , x n , x n , µ n , µ n }. The battery SoC of EV n at time t∆ (t ∈ S) is denoted by σ n (t), σ n (t) :=σ n + ∆ t j=1 x n (j)µ n (j), where µ n (j) = µ n if x n (j) ≥ 0 and µ n (j) = µ n otherwise [2,18]. Accordingly, the SoC profile of EV n is defined by σ n := [σ n (1), . . . , σ n (t), . . . , σ n (T )] ∈ R T ≥0 .

B. EV battery constraints
The EV battery of customer n ∈ N is operated subject to constraints. Specifically, to limit battery over-charging and over-discharging, the SoC is constrained such that σ n ≤ σ n (t) ≤ σ n , for all t ∈ S. The SoC profile, σ n , is then constrained by σ n 1 ≤ σ n ≤ σ n 1, where 1 ∈ R T . To simplify the notation, we define c n := (σ n −σn) ∆ and c n := (σn−σn) ∆ and combine (1) and (2) as c n 1 ≤ Tx n ≤ c n 1, To capture the limited charging/discharging power of the battery, the constraint x n ≤ x n (t) ≤ x n is introduced for all t ∈ S. The charge-discharge profile x n is then constrained by where 1 ∈ R T . To incorporate the EV charging demand e n , where e n := σ * n −σ n (in kWh), the constraint σ n (T ) = σ * n is enforced, which when combined with (1) yields 1 x n = e n /∆, (5) where 1 ∈ R T . To limit the EV charging (and discharging) duration to within the period the EV is grid-connected ([a n , d n ]), a diagonal matrix L n = [l ij ] ∈ R T ×T is defined, in which the i th (1 ≤ i ≤ T ) diagonal entry is 1 if the EV is available across the time-interval (i − 1)∆, i∆ and 0 otherwise. Specifically, if i = j and a n < i ≤ d n , then l ij = 1, otherwise l ij = 0. Clearly, if i = j then l ij = 0. Accordingly, the chargedischarge profile x n is further constrained by where I ∈ R T ×T and 0 ∈ R T . The combined set of battery constraints (3)-(6) for EV n can be succinctly written as Ψ n := {x n : A n x n ≥ b n ,Ā n x n =b n }, where A n := ∈ R T +1 with I ∈ R T ×T , T ∈ R T ×T , 1 ∈ R T and 0 ∈ R T . In (7), the inequality constraint combines (3)-(4) and the equality constraint combines (5)-(6).

C. EV operational costs
A financial policy of net metering, as in [19], is considered, where each customer is billed at the same rate as they are compensated for delivering power to the grid. Daily variations in electricity price are represented by the pricing profile η := [η(1), . . . , η(t), . . . , η(T )] ∈ R T ≥0 , where η(t) is the price of electricity (in $/kWh) over the time period (t − 1)∆, t∆ . Accordingly, the operational cost (in $/day) accrued to customer n ∈ N is defined by where α n is a (very small) regularization constant that yields a smoother battery profile x n with fewer charge-discharge cycles, avoiding unnecessary charging and discharging [20]. The first term in (8) represents the time-varying cost of purchasing (or delivering) energy for EV charging (or discharging). The second term, first introduced in [20], is considered a proxy for the battery degradation cost, which can be alternatively represented by other cost functions as presented in [21]. When Ω n (x n ) > 0, there exists a financial cost to operating the EV, whereas when Ω n (x n ) < 0, there exists a financial benefit.

III. PROBLEM FORMULATION
This section formulates a centralized optimization problem for EV charging and discharging in an unbalanced three phase distribution grid, which serves a set of customers N. Each customer is represented by the REES depicted in Fig. 1.

A. Unbalanced distribution grid model
We consider a three phase, unbalanced, radial distribution feeder with two phase and single phase laterals, represented by the graph G = {V, E}. The set of vertices V := {0, . . . , i, . . . , k} represents nodes along the feeder and the edge-set E ⊆ V × V represents line segments along the feeder, with cardinality |E| = k. Node 0 is the root node (feeder head) that represents an infinite bus, which decouples interactions along the feeder from the wider power grid. The edge (ij) ∈ E exists if there is a line segment between nodes i and j (i, j ∈ V), with node i being closest to the feeder head. Graph G is assumed to be simple with no repeated edges or self loops for any . . , k}, the set of edges on the unique path from node 0 to node i is denoted by Each edge (ij) ∈ E (and each vertex i ∈ V) is either a single phase, two phase, or three phase edge (or node), where phases are labeled a, b and c. The set of phases at node i ∈ V is denoted by Φ i , e.g., Φ i = {a, b, c} for a 3-phase node i. Typically, the infinite bus is represented by a 3-phase node, i.e., Φ 0 = {a, b, c}. The set of phases along edge (ij) ∈ E is denoted by Φ ij , e.g., Φ ij = {a, b, c} for a 3-phase line segment. If edge (ij) exists, then the corresponding phases along edge (ij) are a subset of the phases present at both node i and node j, i.e., Φ ij ⊆ Φ i and Φ ij ⊆ Φ j . In what follows, phases are denoted by φ and ψ, where φ, ψ ∈ {a, b, c}.
Let {(i:φ)} φ∈Φ i represent the tuple of phases at node i ∈ V, e.g., if node i is 3-phase, where |K| = K represents the sum of the total number of phases across all nodes in V. Each element (i:φ) ∈ K represents a supply point that serves N i:φ ≥ 0 number of customers (each represented by the REES in Fig. 1), such that i∈V φ∈Φ i N i:φ = N . We initialize the customer index as n = 1 and we increment the index through all customers at each supply point, and across all supply points in K in ascending order, respectively.
We identify the location of each customer with reference to a supply point, (i:φ) ∈ K. Let Υ denote a binary matrix in which rows are indexed as in K and columns are indexed as in N, such that the element corresponding to row index (i:φ) and column index n is 1 if customer n ∈ N is connected to supply point (i:φ) ∈ K, and 0 otherwise. More formally, where for each (i:φ) ∈ K. Let r ij,φψ and x ij,φψ denote the (self or mutual) resistance and reactance of line segment (ij) ∈ E, respectively, where φ and ψ represent each of the phases along the line segment. Then z ij,φψ := r ij,φψ +ix ij,φψ , so that each line segment has a complex impedance z ij , where for a three-phase line segment Let Z ij,φψ := (îĵ)∈(E i ∩E j ) zîĵ ,φψ ∈ C denote the sum of (self or mutual) impedances along the unique path that is common to the paths from node 0 to node i, and from node 0 to node j, such that both paths incorporate phases φ and ψ. For instance, in Fig. 2 denote a square matrix with row index (i:φ) and column index (j:ψ), where each element in R is defined by R (i:φ), (j: denote a square matrix with row index (i:φ) and column index (j:ψ), where each element in X is defined by X (i:φ), (j: Here, Re{·}, Im{·}, and (·) * denote the real part, the imaginary part, and the conjugate, respectively.
To concatenate voltages across all supply points in K at time t∆, we define Let v 0 denote the squared voltage magnitude of each phase at root node 0, which we set to the squared nominal voltage. Accordingly, let is the net real power injection or consumption arising from Recall, the supply point (i:φ) ∈ K serves N i:φ customers, meaning that all of their REESs with EV loads and non-EV loads connected to supply point (i:φ) contribute to the net real and net reactive power flows, p i:φ (t) and q i:φ (t), respectively. To represent non-EV power flows at a single supply point (i:φ) ∈ K, we denote by p i:φ (t) ∈ R ≥0 (or p i:φ (t) ∈ R <0 ) the net real power (in kW) consumed (or injected) by all non-EV loads connected to phase φ ∈ Φ i at node i ∈ V over the period (t−1)∆, t∆ . Then, we denote by q (i:φ) (t) ∈ R ≥0 (or q (i:φ) (t) ∈ R <0 ) the net reactive power (in kVAR) consumed (or injected) by all non-EV loads connected to phase φ ∈ Φ i at node i ∈ V over the period (t − 1)∆, t∆ . Likewise, to represent EV power flows at a single supply point (i:φ) ∈ K, we denote by ) the net reactive power (in kVAR) consumed (or injected) by all EVs connected to phase φ ∈ Φ i at node i ∈ V over the period (t − 1)∆, t∆ . When real and reactive power losses across the unbalanced distribution grid are ignored as in [7], p i:φ (t) = p i:φ (t) + p i:φ (t) and q i:φ (t) = q i:φ (t) + q i:φ (t) for all t ∈ S. Accordingly, the net real power and net reactive power consumed (and injected) across all supply points in K over the period (t − 1)∆, t∆ are defined by Next, we consider the power flow equations proposed in [23,24], to model the physics of unbalanced, radial distribution grids. Specifically, the authors in [23,24] extend the LinDistFlow equations designed for single-phase grids to the unbalanced case. According to [23,24], the unbalanced version of the LinDistFlow equation is expressed as the squared voltages arising from the aggregate non-EV load, which we refer to as the baseline load at time t∆. Consider the special case where no EV is grid-connected, such that P (t) = 0 and Q(t) = 0. In this special case, the grid voltages (arising from both EV and non-EV loads) reduce to the baseline voltages, i.e., V (t) In what follows, we consider EVs only consume or supply real power (in kW), as done in [2,7]. That is, q i:φ (t) = 0 for all t ∈ S, i ∈ V and φ ∈ Φ i , implying Q(t) = 0. Consequently, we simplify Eq. (11) as, We define by X (t) := [x 1 (t), . . . , x n (t), . . . , x N (t)] ∈ R N the set of charge and discharge rates of all EVs in N over the period (t−1)∆, t∆ , where the n th element represents the charge (or discharge) rate of EV n (n ∈ N). Then, by combining Eq. (9) with X (t), P (t) can be alternatively written as B. Network-aware EV battery scheduling We now formulate a centralized EV charge and discharge optimization problem, in which we constrain the voltage magnitude on each each phase φ ∈ Φ i , at each node i ∈ V, to stay within an operational range [v i:φ , v i:φ ]. Here, v i:φ ∈ R and v i:φ ∈ R denote the lower and upper bound, respectively, for the voltage magnitude on phase φ ∈ Φ i at node i ∈ V.
To constrain the voltages represented by Eq. (14), we To simplify the presentation of constraint (15), we de- (16) In this paper, we seek to minimize the operational cost accrued by each EV customer in N over the planning timehorizon [0, T ∆]. Recall, the operational cost for a single EV customer n, Ω n (x n ), is defined in (8). It follows that we seek to minimize the total operational cost, N n=1 Ω n (x n ), subject to each EV-specific battery constraint set Ψ n in (7) and the simplified voltage constraint in (16). Accordingly, the centralized optimization problem (P1) is succinctly written as, min subject to x n ∈ Ψ n ; ∀n ∈ N,

IV. A DISTRIBUTED APPROACH FOR NETWORK-AWARE EV BATTERY SCHEDULING
In this section, we propose a distributed optimizationbased algorithm to solve the centralized optimization problem (P1), defined in Section III. In the proposed algorithm, EV customers exchange coordination signals via a peer-topeer communication network and determine their EV chargedischarge profiles locally in the absence of a central entity.

A. Communication network
Consider a peer-to-peer communication network represented by an undirected graph G := {N, E}, where the vertex-set N = {1, . . . , n, . . . , N } is the set of EV customers (as defined in Section II) and the edge-set E ⊆ N × N is the set of communication links between EV customers. The presence of edge (nm) ∈ E enables customer n to exchange information with customer m. Accordingly, N n := {m ∈ N | (nm) ∈ E} is defined as the set of neighboring customers from which customer n can receive information. We make the following assumption on graph G. Assumption 1. The undirected graph G := {N, E} is connected, i.e., for every pair of vertices in N, there exists a path of edges in E that connects them.

B. DC-ADMM based network-aware EV battery scheduling
Drawing on the consensus theory [25] and the ADMM formulation [26], we next develop a fully parallelized optimization algorithm that solves (P1) for any connected graph G. Specifically, our development approach is inspired by the DC-ADMM formulation in [17]. The optimization problem (P1) in its present form cannot be solved directly using DC-ADMM. Therefore, we must first transform the inequality constraint (P1.3) into an equality constraint by introducing a vector of slack variables s n ∈ R 2KT for each customer n ∈ N, which satisfy N n=1 Γ n x n + s n = w.
(17) Next, as in [26], we introduce an indicator function I n (u n ) for each customer n ∈ N. Let u n := [x n s n ] ∈ R T +2KT such that u n = [x n (1), . . . , x n (T ), s n (1) . . . , s n (2KT )] , and let S n := {s n : 0 ≤ s n ≤ w}. Then we define I n (u n ) by where Ξ n := Ψ n × S n and Ψ n is as defined in (7). Lemma 1. Let the cost function for a single EV customer n ∈ N be denoted by f n (u n ) := Ω n (u n ) + I n (u n ). Let ξ n := Γ n I ∈ R 2KT ×(T +2KT ) where I ∈ R 2KT ×2KT . Then the optimization problem (P1) is equivalent to the optimization problem (P2), which is defined by min 2) The proof is in the Appendix.
The optimization problem (P2) is a convex minimization problem in which strong duality holds [27]. Therefore, we formulate the dual optimization problem (P3) as follows, min for all n ∈ N. Here, λ ∈ R 2KT is the Lagrange dual variable vector corresponding to constraint (P2.2). The optimization problem (P3) is not separable across customers, since λ is a global vector that needs to be computed by a central entity. Interestingly, by leveraging the DC-ADMM in [17], a distributed counterpart for (P3) can be obtained by equipping each customer n ∈ N with a local estimate (copy) of λ, denoted by λ n ∈ R 2KT . In this way, it is possible to remove the requirement for a central entity and incorporate λ n as a coordination signal amongst EV customers connected via a peer-to-peer communication network. Then, under Assumption 1, (P3) can be equivalently written as (P4): subject to λ n = β nm ; ∀m ∈ N n , n ∈ N, (P4.2) λ m = β nm ; ∀m ∈ N n , n ∈ N, (P4.3) where {β nm } n∈N,m∈Nn ∈ R 2KT is the sequence of consensus auxiliary variable vectors. According to (P4), each customer n ∈ N can minimize its local objective function in (P4.1), subject to consensus constraints (P4.2) and (P4.3).
; ∀m ∈ N n , (19b) where ρ > 0 is a penalty parameter. The initial values of µ nm and γ nm are set as µ According to DC-ADMM, λ [τ ] n is updated by solving the following subproblem for each customer n ∈ N: Then, by substituting Θ n (λ n ) from (P4) into (21), a min-max optimization problem can be obtained, which is given by Since the objective function in (22) is convex in λ n for any u n , and is concave in u n for any λ n , the min-max theorem ( [28], Proposition 2.6.2) is applied to find λ [τ ] n in (22). That is, the subproblem in (22) is solved by its max-min counterpart where the outer maximizer u n and the inner minimizer λ n are determined by Consequently, the solution to subproblem (22) for each customer n ∈ N is obtained by first solving the subproblem (24) with respect to the primal variable u n , followed by evaluating λ n using Eq. (25). Importantly, while DC-ADMM handles the dual problem in (P3), it directly yields the primal optimal solution to (P2), as stated in the following theorem. n } n∈N , as defined in (24), converges to {u n } n∈N .
The proof is in the Appendix.
With the above as background, Algorithm 1 presents our approach to distributed EV charging and discharging, referred to as Dis-Net-EVCD. In Algorithm 1, each of the computational steps is executed by each customer n ∈ N in parallel without requiring the intervention of a central coordinator. Importantly, EVs coordinate amongst themselves by exchanging locally computed λ n vectors, which do not carry EV-specific operational information nor EV-specific battery parameters.
In In more detail, Algorithm 1 consists of a DC-ADMM based iterative routine from lines 8-14 and a receding-horizon routine, as in [2], from lines 17-21. The iterative routine from lines 8-14 is executed until an agreed λ n is reached across all customers in N. We note that the computational steps within the iterative routine (i.e., lines 11-13) require as inputs, quantities that are either locally known by customer n (e.g., Λ n , f n (u n ), ξ n ) or are collected from its neighbors (e.g., {λ   n is updated as per Eq. (25) to proceed to the next iteration. Likewise, the iterative routine continues until a predefined stopping criterion is reached. As proposed in [29], the iterative routine in this paper is set to stop when the iteration number τ exceeds τ max (e.g., τ max = 1000) or = max ν  Inputs: Λn, fn(un) and ξn. 5 Initialize λ Initialize iteration index τ = 0.  n to all neighbors in Nn. 15 until a stopping criterion, as in [29], is satisfied; 16 foreach grid-connected customer n ∈ N do 17 Let un = u  and an = + 1, in preparation for the next time interval. 22 Advance the time index = + 1.
We emphasize that Algorithm 1 is computationally efficient and easy to implement, since line 12 solves a convex optimization problem -for which commercial or open-source software [30] is available -and the updates in lines 11 and 13 involve simple arithmetic computations.

V. NUMERICAL SIMULATIONS
We evaluate Dis-Net-EVCD on the IEEE 13 node test feeder, consistent with the simulation set-up described in [24]. The IEEE 13 node test feeder is available from [22], and it represents an unbalanced distribution circuit with single phase and two phase laterals as depicted in Fig. 2.
We consider charging and discharging of N = 600 customers with charge-only EVs and gridable EVs (V2G-enabled EVs) dispersed across the IEEE 13 node feeder, such that 3, 4, 5, 7, 8, 9, 12}, N Customers are assigned different EV models as described in [31], with battery capacities, c n , ranging from [16,90] kWh. Following [7], each customer n ∈ N is equipped with a level-2 charger supporting a maximum charge rate of x n = 6.6 kW and a maximum discharge rate of x n = −6.6 kW. The remaining design parameters for each EV are as described in [2], where the initial SOC of each EV n ,û n , is within [0.3, 0.5] times the battery capacity, such that there is a uniform distribution across all initial SOCs. For each EV n , σ n = 0.2 × c n , σ n = 0.85 × c n , α n = 5 × 10 −4 , µ n = 0.9 and µ n = 1.1. The EV arrival and departure times, a n and d n , are reflective of survey data [32] gathered by the Victorian Department of Transport in Australia, which contains individual vehicle travel records for a 24 hour period.
The non-EV power consumption, l n , for each customer is obtained from a real-world dataset [33] consisting of residential load collected from households within an Australian distribution grid, as publicly released by Ausgrid, an electricity distribution company in NSW, Australia. The extracted data corresponds to residential energy consumption and solar PV power generation from a set of 300 residential customers on the 5th July 2012, a day on which several peak load conditions were observed. The extracted Ausgrid data is duplicated to populate non-EV power consumption for h = 600 customers. As illustrated in Fig. 4, we consider a ToU financial policy as reported by AGL Energy, an Australian energy provider.
Consistent with [7], a voltage threshold of ±4.6% about the nominal voltage 1 p.u. is considered. All simulations are conducted in Python 3.8.2 + CVXPY [30] on a MacBook Pro with an Intel Core i5 and 8 GB memory.
In Figs. 5(a-b), we present numerical simulations corresponding to a scenario referred to as Price-based EV Charging and Discharging (P-EVCD), where each EV implements a charge-discharge profile designed to minimizing operational costs in the absence of voltage constraints. That is, with P-EVCD, each EV n (n ∈ N) minimizes the operational cost function Ω n (x n ) defined in (8), subject to a local constraint set Ψ n as defined in (7). By contrast, Figs. 5(c-d) present numerical simulations corresponding to Dis-Net-EVCD, where EVs determine charge-discharge profiles as per Algorithm 1, so as to minimize operational costs subject to the voltage constraints and local constraint sets Ψ n of all EVs in N.
In Fig. 5(a) and Fig. 5(c), we present the sum of non-EV power profiles of all customers ( n∈N l n ), referred to as the Aggregate Non-EV Load Profile (ANLP), the sum of EV charge-discharge profiles of all customers ( n∈N x n ), referred to as the Aggregate EV Load Profile (AELP), and the sum of grid power profiles of all customers ( n∈N g n ), referred to as the Aggregate Load Profile (ALP). As illustrated in Fig. 5(a) and Fig. 5(c), the ANLP peaks during the peak-pricing period and then gradually declines through the shoulder-pricing period to the off-peak-pricing period, consistent with the pricing profile η in Fig. 4.
In Fig. 5(b) and Fig. 5(d), we present the voltage profile corresponding to each supply point (i:φ) ∈ K, which is defined by In particular, the grey dotted lines represent the baseline voltage profiles arising from the non-EV load, and black dotted lines represent the grid voltage profiles arising from the combined EV and non-EV load. An important observation from Fig. 5(b) Fig. 4. Pricing profile η taken from [34]: the off-peak period is from 10 pm -7 am, the shoulder-peak period is from 7 am -2 pm and 8 pm -10 pm, and the peak period is from 2 pm -8 pm. Consider P-EVCD as in Figs. 5(a-b). In Fig. 5(a), we observe that the ALP falls to −2.85 MW and rises to 6 MW. The ALP valley corresponding to −2.85 MW occurs at 12 pm when the non-EV load is at its lowest -a consequence of surplus rooftop PV generation. The ALP peak corresponding to 6 MW occurs at 10.30 pm, when EVs collectively charge during the off-peak pricing period. In Fig. 5(a), we also observe that the AELP reflects a short charging period from 12 -2 pm (coincident with shoulder-peak pricing period), arising from gridable EVs (in particular, gridable EVs with lower SoCs). Here, the gridable EV batteries charge ahead of the peak pricing period, in preparation for a discharging cycle during the peak pricing period. We further observe that the AELP continues to rise from 12 -2 pm, as a consequence of more gridable EVs grid-connecting, and thus being available for charging. Then, during the peak period from 2 -8 pm, we observe that the AELP reflects a discharging period, where gridable EVs are compensated to reduce their operational costs. From 2 -8 pm, the AELP gradually descends as a consequence of more gridable EVs grid-connecting, and thus being available for discharging. Then, during the shoulder pricing period from 8 -10 pm, we observe that the AELP reflects a decrease in EV discharge load, since most EVs have an insufficient SoC to continue discharging, and fewer EVs (who have arrived late) continue to discharge.
In Fig. 5(a), we observe that the AELP reflects a substantial charging load during the off-peak pricing period from 10 pm -7 am, which fills the over-night valley of the ANLP. As a result, in Fig. 5(b), we observe that the grid voltage profiles arising from P-EVCD (black dotted curves) undergo voltage dips that exceed the lower voltage boundary -resulting in voltage violations across the IEEE 13 node distribution feeder. Specifically, the largest voltage violation in Fig. 5(b) occurs at 10.30 pm, coincident with the highest ALP peak in Fig. 5(a), which is caused by the enormous charging load during the off-peak pricing period.
Consider Dis-Net-EVCD as in Figs. 5(c-d). In Fig. 5(c), we observe that the ALP corresponding to Dis-Net-EVCD is flatter (with less prominent peaks and valleys) when compared to the ALP from P-EVCD in Fig. 5(a). Specifically, the ALP in Fig. 5(c) has a peak of 3.66 MW at 10.30 pm, which is 39% less than the respective ALP peak in Fig. 5(a). Such a peak load reduction is a result of the voltage constraint incorporated in Dis-Net-EVCD, which restricts EV charging and discharging to levels that maintain voltages within upper and lower limits. In Fig. 5(d), we observe that when EVs implement Dis-Net-EVCD, the grid voltage profiles (black dotted curves) across all supply points in the feeder remain within the upper and lower limit of ±4.6% about the nominal voltage. Importantly, Dis-Net-EVCD, as opposed to P-EVCD, corrects the voltage violations in baseline voltage profiles as observed from 12 -12.30 pm and from 5 -9.30 pm, and regulates the grid voltage profiles to be within the operating voltage envelope across the day. In addition, Dis-Net-EVCD completes all EV charging demands ahead of their departures.
The total operational cost, N n=1 Ω n (x n ), corresponding to P-EVCD in Figs. 5(a-b) is 391 AUD, and Dis-Net-EVCD in Figs. 5(c-d) is 483 AUD, meaning that voltage regulation with Dis-Net-EVCD results in an increased operational cost compared to P-EVCD. However, Dis-Net-EVCD significantly reduces operational costs (by 78%) when compared to uncoordinated EV charging (where EVs are assumed to be charging at their full power immediately upon arrival), in which case the operational cost is 2173 AUD.
The computational time required to produce the results depicted in Figs. 5(c-d) was approximately 5.4 s for each time index of the receding time horizon. The respective computational time excluded latencies associated with peerto-peer communication. A centralized CVXPY Python solver, would otherwise take 5.5 minutes to produce the same results. That is, Dis-Net-EVCD is approximately 60 times more computationally efficient when compared to its centralized counterpart, where a central entity is assumed to be computing the charge-discharge profiles for all EVs.
In Fig. 6, we compare the convergence of Dis-Net-EVCD with an alternative distributed and iterative algorithm called Distributed Lagrangian Primal-Dual Subgradient (DLPDS), proposed in [35]. In Fig. 6, the error is the normalized objective function gap between the centralized (global) algorithmic solver and the distributed algorithmic solver, corresponding to a respective iteration. Formally, error = (obj−obj * ) obj * , where obj * denotes the objective function (P1.1) result obtained with the centralized solver and obj denotes the objective function (P1.1) result obtained with the distributed algorithm at a certain iteration. With Dis-Net-EVCD, we observe that the error reduces to 10 −5 within only 45 iterations. In fact, the speed of convergence with Dis-Net-EVCD is much faster than DLPDS. Moreover, DLPDS also requires the communication network to be strongly connected and doubly stochastic [35], whereas Dis-Net-EVCD simply requires the communication network to be connected.
In Fig. 7, we quantify the error introduced through linearization of the unbalanced IEEE 13 node feeder, by comparing the grid voltage profiles computed with the non-linear DistFlow equations derived in [36], against the LinDistFlow equations derived in [23,24]. Specifically, we apply the EV chargedischarge profiles returned by Dis-Net-EVCD together with the residential load profiles, l n , to compute the grid voltage profiles for the linearized and non-linearized representations of IEEE 13 node feeder. In Fig. 7, we observe that the voltage gap between the linearized and non-linearized representations is as low as 0.004 p.u. Importantly, by including a lower voltage constraint of 0.954 p.u. within Dis-Net-EVCD, we observe that the voltage drop across the IEEE 13 node feeder is limited to 0.95 p.u. with a DistFlow representation of the power flow equations. A steady-state voltage of 0.95 p.u is within the voltage operating standard, as reported in [7,37]. That is, by Fig. 6. Convergence of Dis-Net-EVCD against DLPDS [35]. The error versus iteration number. Fig. 7. Black and grey dotted lines are grid voltage profiles obtained using the LinDistFlow model [23,24] and the DistFlow model [36], repectively.
reducing the voltage operating envelope (i.e., by reserving a buffer voltage zone) within Dis-Net-EVCD, it is possible to correct for voltage deviations corresponding to line losses that are otherwise omitted in the LinDistFlow equations.

VI. CONCLUSION
We have presented a distributed optimization-based, network-aware algorithm called Dis-Net-EVCD to coordinate EV charging and discharging in unbalanced distribution grids. Specifically, Dis-Net-EVCD reduces the operational costs incurred by EV customers, while regulating voltages across the distribution grid to be within the operational limits. We have proved that Dis-Net-EVCD converges to the global optimum of the underlying optimization. Through numerical simulations, it was demonstrated that Dis-Net-EVCD achieved a peak load reduction of 39% compared to a benchmark method of EV charging that minimized customer operational costs in the absence of voltage constraints. Importantly, Dis-Net-EVCD satisfied all EV charge requirements ahead of their intended departure times while alleviating potential under-and over-voltage conditions, and thereby improving power quality across the distribution grid. Future work to extend Dis-Net-EVCD to incorporate power flow limits and dynamic thermal grid constraints, is possible.

A. Proof of Lemma 1
The result directly follows by substituting voltage constraint (17) and indicator function (18) in the optimization problem (P1).
To show (27), we consider the optimality condition of (24),  (33) verifies that the statement (28) is also true. This completes the second part of the proof.