Zoning Evaluation for Voltage Optimization in Distribution Networks with Distributed Energy Resources

This paper deals with the problem of the voltage profile optimization in a distribution system including distributed energy resources. Adopting a centralized approach, the voltage optimization is a non-linear programming problem with large number of variables requiring a continuous remote monitoring and data transmission from/to loads and distributed energy resources. In this study, a recently-proposed Jacobian-based linear method is used to model the steady-state operation of the distribution network and to divide the network into voltage control zones so as to reformulate the centralized optimization as a quadratic programming of reduced dimension. New clustering methods for the voltage control zone definition are proposed to consider the dependence of the nodal voltages on both active and reactive powers. Zoning methodologies are firstly tested on a 24-nodes low voltage network and, then, applied to the voltage optimization problem with the aim of analyzing the impact of the R/X ratios on the zone evaluation and on the voltage optimization solution.


Introduction
In the normal operation of modern distribution systems, variations of active and reactive powers absorbed and/or injected by controllable loads, distributed generators (DGs), and energy storage systems (ESSs) as well as changes in network topology can cause violation of the nodal voltage limits, fixed to ±10% of the rated voltage of the grid.
Presently, the utility-side voltage regulation cannot adequately respond to the voltage variations. In fact, the conventional Volt/Var control devices (i.e., on-load tap changer of the transformer in the HV/MV substations, step voltage regulators and capacitor banks) have slow response time and can negatively interfere with distributed energy resources (DERs) (i.e., DGs, ESSs including electrical vehicles, and smart loads), deteriorating voltage regulation. To keep the nodal voltages within the admissible limits, new flexible voltage control strategies have to be developed, involving both conventional devices and DERs [1].
The recently-proposed architectures for voltage regulation can be classified into local control and communication-based control according to whether the use of communication infrastructures is not or is required. The second category can be further split into: centralized, distributed and decentralized control architectures [2].
In a local control, the controllers of both DERs and conventional Volt/Var devices use measurements at the point of common coupling (PCC) to individually elaborate and actuate a control action without requiring any information exchange. The result is a globally non-optimal control solution. Undoubtedly, this is the lowest-cost solution to start involving DERs in the voltage control of existing distribution networks. Unfortunately, the absence of coordination may introduce technical problems related to the interaction among controllers or even system instability [3,4].
In a centralized control, a central control unit, located at the substation level, firstly solves a system voltage optimization/control problem by receiving measurements from all the nodes of the network; then, it sends the optimal set-points back to the local controller of both DERs and conventional Volt/Var devices. This approach yields the optimal solution but its implementation is very expensive requiring a large communication infrastructure with adequate bandwidth to exchange information quickly and accurately [5][6][7][8].
In a distributed control, the controllers of both DERs and conventional Volt/Var devices use measurements at PCC to achieve local voltage regulation as in a local control but they also exchange information with the controllers of the neighboring nodes. This approach aims at overcoming the technical problems of the local control while limiting the investments for the communication infrastructure [9,10].
In a decentralized approach, the distribution system is usually partitioned into voltage control zones (VCZs). Each VCZ contains nodes weakly coupled with the nodes belonging to other VCZs from an electrical point of view; in the 'electric center' of a VCZ is placed the pilot node (PN), whose voltage variation best represents the variation of the voltage in the VCZ. In each VCZ a centralized control is present whereas the coordination is obtained by a distributed control among the PNs. This approach guarantees the optimal solution for each VCZ and it requires a communication infrastructure lighter than the centralized control [11,12].
The division of the network in VCZs require the definition of an electrical distance, a clustering criterion and, in some case, a clustering quality index [13]. In literature, various measures of the electrical distance have been proposed, based on geographical distance, line impedance, voltage sensitivity coefficients [14] and the most commonly used involves the sensitivity of the nodal voltages to the reactive power injections [13,[15][16][17][18][19][20][21]. Such an approach is typically adopted for the secondary voltage control in transmission systems. Unfortunately, in distribution systems nodal voltages depends on both active and reactive power injections due to the high R/X ratio of the grid. To face this problem, in [22] electrical distances are evaluated taking into account the dependence of the nodal voltages to the active powers. In [11,12] both the sensitivities of the nodal voltages to active and reactive powers injected by DERs are considered and used for the definition of a weighted linear electrical distance. In [23] a modified electrical distance is proposed combining the voltage sensitivities in [11,12] with the conventional definition of electrical distances derived from the line impedance measure. Moreover, there are also many methods for grouping the nodes with similar electrical distances in VCZs, such as hierarchical clustering [11,13,16,17], K-means clustering [18], genetic algorithms [19,24].
In this paper, extending and improving the approach proposed in [14], a partition of the distribution network in VCZs that account for the sensitivities of the nodal voltages to active and reactive powers is firstly proposed and then applied to optimize the voltage profiles of the distribution feeders with DERs. In particular, using the sensitivity coefficients derived from the linear method for the steady-state analysis of the distribution system in [25,26], various definitions of the electrical distance are proposed in this paper, which consider separately or simultaneously the dependence of the nodal voltages on active and reactive power injections, independently from the number and position of the DERs. Successively, the hierarchical clustering algorithm in [16] is adopted and extended to adequately work also with two electrical distances defined for each node. Then, adopting the centralized approach as in [15], the zoning methodologies are applied to optimize the voltage profile of the distribution systems with DERs so as to guarantee a solution close to the optimal one and to reduce the requirements for the monitoring and communication infrastructures. Finally, the influence of different R/X ratios on the proposed zoning methodologies as well as on the voltage optimization problem is analyzed in a case study. Differently from [11,12,22,23], in this paper various electrical distance measures have been proposed accounting for active and reactive power sensitivities, which are identified on the basis of the structural characteristics of the power grid thanks to the linear model in [25,26]. Compared to [15], in this paper: (i.) the sensitivities of the nodal voltages to not only reactive but also active powers are considered; (ii.) the hierarchical clustering algorithm is extended to obtain aggregation based on active and reactive powers; (iii.) the centralized voltage optimization problem includes DERs with controllable active and reactive powers; (iv.) the impact of the R/X ratios on the VCZ evaluation and on the voltage optimization solution is analyzed and general considerations are derived.
The paper is organized as follows: the linear method for the steady-state operation of the distribution system is briefly recalled in Section 2. Section 3 presents the proposed zoning methodologies and includes Appendices A and B containing the extensions of the hierarchical clustering algorithms. In Section 4 the application of the zoning methodology to the voltage optimization problem is reported. Referring to a 24-nodes distribution system, in Section 5 various network partitioning options are tested and, then, applied to the centralized voltage optimization problem.

Brief Recalls of the Linear Method
A Jacobian-based method for the steady-state analysis of radially-operated distribution systems with DERs has recently been presented in [25,26]. Starting from an initial operating point, it provides the closed-form analytical expressions of the sensitivity coefficients that linearly relate the variations of the electrical variables characterizing each node of a feeder to the variations of the powers injected by all the DERs connected to the grid.
Referring all the variations with respect to the initial operating point, the steady-state linear model of a radial distribution network with N nodes including DERs is expressed as: where: • ∆P i and ∆Q i are the variations of, respectively, the active and the reactive powers out-flowing the i-th node of the grid; • ∆V 2 i is the variation of the squared nodal voltage at the i-th node of the grid; • K der is a subset of {1, . . . , N} including only the nodes of the network with DERs; • ∆P der k and ∆Q der k are the variations of, respectively, the active and reactive powers injected by a DER at the k-th node of the grid.
The sensitivity coefficients in (1), relating the electrical variables at the i-th node to the power injections of the DER connected at the k-th node, can be collected in the (3 × 2) matrices T i (k): so that: In Equations (1) and (3) DERs are assumed to be equipped with P-Q control; the extension of the model to the distribution network including DERs with different types of controls is detailed in [15]. It is worth noting that the adoption of the new variable ∆V 2 i (in place of ∆V i ) permits to limit the approximation caused by the linearization of the power flow equations with respect to other Jacobian-based methods.

Zoning Methodologies
The aim of a zoning methodology is to obtain a simplified representation of the distribution network suitable for voltage control. Basically, the problem involves four steps [13]:

•
Defining an electrical distance among the nodes of the network; • Grouping the nodes into VCZs; • Identifying a PN in each VCZ; • Determining the optimum number of VCZs.
The electrical distance among the nodes as well as the topology of the network are important to deeply investigate the structure of a power grid. The electrical distance is a measure of the electrical connectedness among the nodes. Various measures of the electrical distance have been proposed [14] but the most commonly used involves the Jacobian matrix of the load flow equations and its inverse, named sensitivity matrix. Adopting and extending the approach outlined in [15], various definitions of the electrical distance are proposed in this paper, based on the sensitivity coefficients derived from the linear model recalled in Section 2 with respect to both active and reactive powers.
After the first step, nodes with similar electrical distances are grouped to form a VCZ; those nodes in turn result to be weakly coupled from an electrical point of view with the nodes belonging to other VCZs. In other words, the problem of grouping the nodes in VCZs is transformed into a clustering problem in the space of the electrical distance. Clustering can be achieved by various algorithms [20] but one of the most commonly used is the hierarchical algorithm [11,13,16,17], which is adopted in this paper and extended to account for various electrical distances.
Once the VCZs are available from the second step, for each VCZ a PN is identified. The PN is chosen as the node whose voltage variation best represents the voltage variations of all the nodes in the VCZ. The identification of the PN of a VCZ is based on an extension of the algorithm in [20] that evaluates the proximity of each node to all the other nodes belonging to the same VCZ.
Finally, in the last step, the best number of VCZs is chosen. To this aim, an index that represents the goodness of the clustering results must be used; in this paper, the Silhouette index [27,28] is used because it is the most adequate for hierarchical algorithms.

Electrical Distances
In a radial distribution network with a DER connected at the k-th node, the square nodal voltage at the i-th node varies for the changes of both DER active and reactive power injections according to: The voltage sensitivity coefficients to the DER active and reactive powers can be evaluated by extracting the third row of T i (k) in (2): Assuming to connect a DER in each node of the network (i.e., k = 1, . . . , N), the voltage sensitivity coefficients to the DER active and reactive powers can be collected in two (N × N) voltage sensitivity matrices, respectively Γ P and Γ Q , whose generic elements are: The electrical distances among the nodes of the network can be evaluated by using three sensitivity measures described in the following, namely: • Sensitivity to the only active power; • Sensitivity to the only reactive power; • Sensitivity to both active and reactive powers.

Sensitivities to Active Powers
Considering only the impact of the active power, the variation of the square nodal voltage at the i-th node to the change of the DER active power at the k-th node is equal to: Similarly, the variation of the square nodal voltage at the k-th node to the change of the DER active power at the k-th node is equal to: Dividing Equation (9) by Equation (10), the variation of the voltage at the i-th node to the change of the DER active power at the k-th node can be rewritten in terms of 'voltage attenuation' of the voltage variation between the nodes i and k caused by an active power variation in k as: The concept of electrical distance used in this paper is presented in [16]. To define a measurement with symmetrical property, the product of the voltage attenuations between nodes i and k and between nodes k and i is considered. Moreover, to obtain a null distance between a node and itself the logarithmic function is used. The electrical distance between the nodes i and k, D P i,k , is then defined as: The values of the electrical distances are usually normalized according to: and, then, collected in the (N × N) electrical distance matrix D P ; such a matrix represents a network divided in N VCZs.

Sensitivities to the Reactive Powers
Considering only the impact of the reactive power, the electrical distance between the nodes i and k, D Q i,k is defined in a similar way as for the active powers, see Equation (12): and the normalized measures are collected in the (N × N) electrical distance matrix D Q .

Sensitivities to Active and Reactive Powers
Considering the impact of both active and reactive powers on the voltage amplitudes, the electrical distance between the nodes i and k can be defined using the mathematical concept of Manhattan distance as: or adopting the Euclidean distance: After the normalization D 1 are collected into the (N × N) electrical distance matrices D 1 PQ and D 2 PQ , respectively.

Voltage Control Zones
The aggregation of the nodes in VCZs is performed by using the previously-defined electrical distances. The related matrices (D P , D Q , D 1 PQ and D 2 PQ ) represent a starting point, that is a network composed of N VCZs, each VCZ coinciding with a node. The aim is to aggregate the N nodes into VCZs yielding the set VCZ = {VCZ 1 , . . . , VCZ N vcz } composed by the N vcz subsets VCZ h (with h = 1, . . . , N vcz ) containing the nodes of the network grouped in each zone. The clustering is carried by applying a hierarchical algorithm so as to perform: • Aggregation based on the only active power; • Aggregation based on the only reactive power; • Aggregation based on both active and reactive powers.

Aggregation Based on Active Power
The electrical distances take into account the dependence of the voltage amplitudes on the only active power injections. Then, reference is made to the matrix D P . The main steps of the zoning algorithm, hereafter referred to as P algorithm, are summarized in [15]. The outputs of the zoning algorithm are the N vcz subsets VCZ P h with h = 1, . . . , N vcz containing the nodes of the network grouped in each VCZ.

Aggregation Based on Reactive Power
The electrical distances take into account the dependence of the voltage amplitudes on the only reactive power injections with the matrix D Q . The algorithm, hereafter referred to as Q algorithm, is the same as in the previous case except for the distance matrix and the outputs are the subsets VCZ Q h with h = 1, . . . , N vcz containing the nodes grouped in each VCZ.

Aggregation Based on Active and Reactive Powers
The electrical distances take into account the dependence of the voltage amplitudes on both active and reactive power injections. Four algorithms are considered.
In the first algorithm, hereafter referred to as PQ algorithm, the hierarchical algorithm in [15] is extended to coordinately work on both matrices D P and D Q . The main steps of the algorithm are detailed in Appendix A. The output of the zoning algorithm are the subsets VCZ PQ h with h = 1, . . . , N vcz .
In the second algorithm, hereafter referred to as P∩Q algorithm, the VCZs are determined as the intersection sets of the N vcz VCZs obtained by separately applying the two hierarchical algorithms in Sections 3.2.1 and 3.2.2. The main steps of the intersection algorithm are detailed in Appendix B. The output is composed of N vcz subsets VCZ P∩Q h with h = 1, . . . , N vcz , being N vcz ≥ N vcz .
In the third and fourth algorithms (hereafter referred to as D 1 and D 2 algorithms, respectively), the hierarchical algorithm described in [15] is applied to the matrices D 1 PQ and D 2 PQ , respectively, providing as output the subsets VCZ PQ 1 h and VCZ PQ 2 h with h = 1, . . . , N vcz , respectively.

Pilot Nodes
The evaluation of the PN in each VCZ is performed by applying the same algorithm whichever the chosen type of aggregation, with the only exception of the PQ and the P ∩ Q algorithms. Referring to the generic matrix D of electrical distances and to the subset of nodes clustered in the generic h-th zone VCZ h (with h = 1, . . . , N vcz ) and composed of m h nodes, the PN is determined according to the following steps: For the PQ and the P∩Q algorithms, the submatrices D h vcz and the minimum sum must be evaluated for both the matrices D P and D Q . The PN is chosen as the node corresponding to the row with the minimum sum, among the rows of both the matrices.
In the following, the set of the N vcz PNs is indicated as H pn , which is the subset of {1, . . . , N} including the indices of the nodes that are PNs.

Clustering Quality Indices
To evaluate the quality of the performed clustering, it is useful to adopt a clustering quality index. A first distinction among clustering indices can be done based on the type of the clustering algorithm, which is either hierarchical or partitional [24,25]. The algorithm used in this paper belongs to the first category, namely an agglomerative hierarchical clustering. The most used index for hierarchical clustering is the Silhouette index.
Let the set of the indices of the nodes belonging to VCZ h be represented as n h 1 , . . . , n h m h and let D vcz (n h i , n h k ) represent the distance between the nodes n h i and n h k . The average distance D vcz (n h i ) between the node n h i and the other nodes in the same cluster VCZ h is given by: whereas the average distance D k vcz (n h i ) among the node n h i and all the other nodes in another clusters VCZ k (with k = 1, . . . , N vcz and k = h) is given by: The minimum average distance D min vcz (n h i ) among the node n h i and all the other nodes in the other clusters different from VCZ h is given by: The Silhouette index s(n h i ) of the node n h i in the cluster VCZ h is defined as: and bounded within the range [−1, 1]. It is possible to define the Silhouette index S h vcz of the VCZ h as as well as of the whole clustering S as: The ; that means that the average distance from nodes of the same cluster is larger than the one from the nodes of at least one other cluster. In general, the closer the Silhouette index is to 1, the better the clustering has been performed.

Application to a Voltage Optimization Problem
The voltage control is structured in a hierarchical architecture that incorporates three control levels [20]. The voltage optimization problem is implemented at the secondary control level. In a centralized approach, the voltage optimization is typically implemented on a specific platform in the substation, called Distribution Management System (DMS). It is formulated as a minimization problem subject to the power flow (PF) equations and network constraints. The solution of the optimization is periodically provided (typically with a time interval of some minutes) and determines the set-points of the controllers of the conventional Volt/Var control devices as well as of the DERs connected to the grid [8,29,30].
In this paper, the optimization problem minimize the sum of the squared distances of the, squared nodal voltages from their reference values acting only on the active and reactive power provided by DERs according to: where: P der min,k , P der max,k and Q der min,k , Q der max,k are the available ranges of, respectively, the active and reactive powers of the k-th DER.
The voltage control is achieved by:

•
Measuring the voltage amplitude at the slack node; • Collecting the measures of active and reactive powers injected/absorbed by DERs and loads; • Solving a minimization problem subject to PF equations and network constraints; • Sending the optimal set-points to the local controllers of the DERs installed in the grid; as highlighted in Figure 1a. Equation (23) is a non-linear programming problem of large dimension. Despite it provides the best solution of the voltage control problem, such an approach presents the following two drawbacks. Firstly, possible convergence problems may arise in the numerical solution, due to the non linear PF equations and to the large dimension; secondly, expensive communication infrastructure is needed to collect real-time measures from all DERs and loads. the set-points of the controllers of the conventional Volt/Var control devices as well as of the DERs connected to the grid [8,29,30]. In this paper, the optimization problem minimize the sum of the squared distances of the, squared nodal voltages from their reference values acting only on the active and reactive power provided by DERs according to: where: is the acceptable variation of ; are the available ranges of, respectively, the active and reactive powers of the k-th DER.
The voltage control is achieved by:

•
Measuring the voltage amplitude at the slack node; • Collecting the measures of active and reactive powers injected/absorbed by DERs and loads; • Solving a minimization problem subject to PF equations and network constraints; • Sending the optimal set-points to the local controllers of the DERs installed in the grid; as highlighted in Figure 1a. Equation (23) is a non-linear programming problem of large dimension.
Despite it provides the best solution of the voltage control problem, such an approach presents the following two drawbacks. Firstly, possible convergence problems may arise in the numerical solution, due to the non linear PF equations and to the large dimension; secondly, expensive communication infrastructure is needed to collect real-time measures from all DERs and loads.
(a) (b) To reduce the number of real-time measurements and to avoid using non linear PF equations, the linear modeling of the distribution system and its clustering into VCZs with PNs can be exploited as in [15]. Then, the minimization problem is rewritten in a variation form with respect to the initial operating point as: To reduce the number of real-time measurements and to avoid using non linear PF equations, the linear modeling of the distribution system and its clustering into VCZs with PNs can be exploited as in [15]. Then, the minimization problem is rewritten in a variation form with respect to the initial operating point as: where: • In this way, the centralized voltage control can be achieved by: • Measuring the voltage amplitude at the slack node; • Collecting only the measures of the voltage amplitudes at the PNs; • Solving a minimization problem subject to linear model (The linear model in (3) is reduced to the variations of the nodal voltage amplitudes) and network constraints; • Sending the optimal set-points to the local controller of the DERs installed in the grid as highlighted in Figure 1b.
The main positive features of this formulation are: (i.) it does not present convergence problem due to the linearity of the model of the distribution system; (ii.) it requires a limited communication system because only the voltage amplitudes at the PNs are collected; (iii.) the convergence of the solving algorithm is guaranteed and the solution requires low computational burden, because it is a quadratic programming problem of small dimension. As a drawback, the solution of Equation (24) is an approximation of the solution of Equation (23) but very close to the optimal one, as shown in the case study, provided that the network has adequately been clustered.

Results
Reference is made to the 24-nodes LV distribution network in Figure 2a. A 20/0.4 kV substation with a 25 kVA transformer feeds a network with 22 nodes composed of 2 main feeders, each one with 2 laterals. Electrical parameters of lines and uncontrolled loads are reported in Tables 1 and 2, respectively. In the following the values in p.u. are referred to a power basis equal to 25 kVA. Results are organized as follows. Zoning methodologies are firstly tested in the 24 nodes LV network with the aim to highlight their specific characteristics. Then, the applicability of the zoning techniques to the voltage optimization problem is investigated, with the aim to obtain a simplified representation of the LV network suitable for smart control of six DERs connected to the network, see Figure 2b.

Zoning Methodologies
The linear method is applied to evaluate the full (3 × 2) sensitivity matrices ( ) with , = 3, … , 24 in the following assumptions: • No-load voltage amplitude at the MV busbar: = 1.0 p.u.; • Active and reactive powers absorbed by uncontrolled loads equal to 70% of the rated values: = 70% , and = 70% , ; • No DER connected to the grid.
Then, the sensitivity coefficients of the square nodal voltages to the DER active and reactive powers injections are extracted from the (3 × 2) sensitivity matrices ( ) with = 3, … , 24 (third row) and collected in the (22 × 22) voltage sensitivity matrices, respectively, and , once and for all. The (22 × 22 ) matrices (0) of the electrical distances between any two nodes of the network are evaluated by the sensitivity measures reported in the following: • Sensitivity to active power, referred to as (0); • Sensitivity to reactive power, referred to as (0); • Sensitivity to active and reactive powers based on the Manhattan distance, referred to as (0); • Sensitivity to active and reactive powers based on the Euclidean distance, referred as (0).
Nodes are merged in VCZs by applying the six different methodologies listed below (see Section 3): • Hierarchical clustering with (0), referred to as P algorithm; • Hierarchical clustering with (0), referred to as Q algorithm; • Hierarchical clustering with (0) and (0), referred to as PQ algorithm; • Intersection clustering with (0) and (0), referred to as P ∩ Q algorithm; • Hierarchical clustering with (0), referred to as D 1 algorithm; • Hierarchical clustering with (0), referred to as D 2 algorithm. and PNs are derived by using the algorithm reported in Section 3.3. Assigned the topology of the distribution network, for each zoning methodology, the three cases described in the following are considered to analyze the dependence of the results on the R/X ratio of the line parameters: • Case R-In this base case, the R/X ratio of the line parameters are reported in Table 1. They are typical values of a LV network in which resistive parameters are prevailing; note that the different

Zoning Methodologies
The linear method is applied to evaluate the full (3 × 2) sensitivity matrices T i (k) with i, k = 3, . . . , 24 in the following assumptions: Then, the sensitivity coefficients of the square nodal voltages to the DER active and reactive powers injections are extracted from the (3 × 2) sensitivity matrices T i (k) with i = 3, . . . , 24 (third row) and collected in the (22 × 22) voltage sensitivity matrices, respectively, Γ P and Γ Q , once and for all. The (22 × 22) matrices D(0) of the electrical distances between any two nodes of the network are evaluated by the sensitivity measures reported in the following: • Sensitivity to active power, referred to as D P (0); • Sensitivity to reactive power, referred to as D Q (0); • Sensitivity to active and reactive powers based on the Manhattan distance, referred to as D 1 PQ (0); • Sensitivity to active and reactive powers based on the Euclidean distance, referred as D 2 PQ (0). Nodes are merged in VCZs by applying the six different methodologies listed below (see Section 3):

•
Hierarchical clustering with D P (0), referred to as P algorithm; • Hierarchical clustering with D Q (0), referred to as Q algorithm; • Hierarchical clustering with D P (0) and D Q (0), referred to as PQ algorithm; • Intersection clustering with D P (0) and D Q (0), referred to as P ∩ Q algorithm; • Hierarchical clustering with D 1 PQ (0), referred to as D 1 algorithm; • Hierarchical clustering with D 2 PQ (0), referred to as D 2 algorithm. and PNs are derived by using the algorithm reported in Section 3.3. Assigned the topology of the distribution network, for each zoning methodology, the three cases described in the following are considered to analyze the dependence of the results on the R/X ratio of the line parameters:

•
Case R-In this base case, the R/X ratio of the line parameters are reported in Table 1. They are typical values of a LV network in which resistive parameters are prevailing; note that the different ratios are related to the presence of both lines and cables (R/X is higher for cables) and to difference sections (R/X is lower at the head of the feeder).

•
Case R≈X-Each R/X ratio in Table 1 is multiplied for a factor equal to 0.2 so as to obtain comparable values of resistive and inductive parameters. • Case X-Each R/X ratio in Table 1 is multiplied for a factor equal to 0.05; results are typical of a MV network in which inductive parameters are prevailing.
The results obtained by clustering the nodes of the distribution network into 3 VCZs are shown in Table 3, reporting the nodes assigned to each VCZ and the corresponding PN (which are obtained by using the six different clustering algorithms) for the three considered cases of R/X ratios. In the Case R, most of the algorithms taking simultaneously into account the sensitivities of the nodal voltages to the active and reactive powers (i.e., PQ, D 1 and D 2 ) gives rise to the same VCZs and PNs derived from the application of the P algorithm; this results was expected, since the network is predominantly resistive and, then, more sensitive to the active power injections. In the Case R≈X, VCZs and PNs are identical for all the clustering algorithms. In the Case X the result is opposite to the one in the Case R: since the network is predominantly inductive, all the algorithms that aggregate the sensibilities to the active and reactive powers (i.e., PQ, D 1 and D 2 ) provide the same results of the Q algorithm; it is worth noticing that in this case the P∩Q algorithm does not provide an intersection of the network into three VCZs.  Tables 4-7 for a network divided, respectively, in 4-7 VCZs. Referring to Table 4, both VCZs and PNs are the same for P, D 1 and D 2 algorithms as well as for Q and PQ algorithms, whatever the value of the R/X ratio; the P∩Q algorithm does not provide a solution in the Case R≈X. Although the clustered nodes of the VCZs are different, the resulting PNs are quite similar among the six algorithms and the three cases. The outcomes of the clustering of the network into five VCZs, reported in Table 5, exemplifies the most classical results. In the Case R, that is a network with dominant resistive parameters, all the algorithms follow the same clustering provided by the P algorithm, with the only exception of the Q algorithm; in the dual Case X, that is a network with dominant resistive parameters, the result is the opposite and all the algorithms follow the same clustering provided by the Q algorithm, with the only exception of the P algorithm; eventually, in the intermediate Case R≈X all six algorithms provide the same VCZs and PNs. Table 6 shows the VCZs and PNs in the case of a network clustered into six VCZs; surprisingly, results are the same for all the six algorithms and all the three cases. For the sake of completeness, Table 7 reports the results obtained by clustering the network into seven VCZs; indeed, they are of little significance because there are VCZs with only one node, due to the large number of VCZ with respect to the total number of network nodes.
With reference to the P and the Q algorithms, the Silhouette index was calculated and plotted with respect to the number of VCZs in Figure 3, for the three considered cases. Comparing the first two plots in Figure 3 which refer to the Case R, the index is generally higher for the P algorithm as expected, because R is prevalent with respect to X; for the P algorithm the index increases in value with the increasing number of VCZ and the highest values are reached for VCZ ≥ 5; for the Q algorithm the index has a slight peak for VCZ equal to 4. The second couple of plots refer to the Case R≈X; the index for the Q algorithm starts from a lower value for VCZ equal to 2 with respect to the previous Case R, then it increases rapidly to reach a peak for VCZ equal to 4 as previously; the index for the P algorithm is a bit smaller than the one for the Case R, but starting from a similar value for VCZ equal to 2, reaches its peak for VCZ equal to 5, then it decreases. In the last Case X, the index for the Q algorithm has a similar behavior as the one for the P algorithm in the first Case R, that was expectable; in fact, the index increases with the number of VCZ to get the best values for VCZ ≥ 5. The index for the P algorithm has a peak for VCZ equal to 3, then it decreases with increasing number of VCZs. Some conclusions can be drawn:

•
In the Case R, index for the P algorithm shows higher values and the best choice is for VCZ ≥ 5; • In the Case R≈X, both P and Q algorithms show the largest index value for VCZ = 5, considering that for the Q algorithm there is small variation of the index between VCZ equal to 4 or 5; • In the Case X, the index for the Q algorithm suggests a clustering composed of VCZ equal to 5 or larger values whereas the index for the P algorithm has a peak for VCZ equal to 3.
With reference to the P and the Q algorithms, the Silhouette index was calculated and plotted with respect to the number of VCZs in Figure 3, for the three considered cases. Comparing the first two plots in Figure 3 which refer to the Case R, the index is generally higher for the P algorithm as expected, because R is prevalent with respect to X; for the P algorithm the index increases in value with the increasing number of VCZ and the highest values are reached for VCZ ≥ 5; for the Q algorithm the index has a slight peak for VCZ equal to 4. The second couple of plots refer to the Case R≈X; the index for the Q algorithm starts from a lower value for VCZ equal to 2 with respect to the previous Case R, then it increases rapidly to reach a peak for VCZ equal to 4 as previously; the index for the P algorithm is a bit smaller than the one for the Case R, but starting from a similar value for VCZ equal to 2, reaches its peak for VCZ equal to 5, then it decreases. In the last Case X, the index for the Q algorithm has a similar behavior as the one for the P algorithm in the first Case R, that was expectable; in fact, the index increases with the number of VCZ to get the best values for VCZ ≥ 5. The index for the P algorithm has a peak for VCZ equal to 3, then it decreases with increasing number of VCZs. Some conclusions can be drawn: • In the Case R, index for the P algorithm shows higher values and the best choice is for VCZ ≥ 5; • In the Case R≈X, both P and Q algorithms show the largest index value for VCZ = 5, considering that for the Q algorithm there is small variation of the index between VCZ equal to 4 or 5; • In the Case X, the index for the Q algorithm suggests a clustering composed of VCZ equal to 5 or larger values whereas the index for the P algorithm has a peak for VCZ equal to 3.

Voltage Optimization Problem
To test the effectiveness of dividing the distribution network into VCZs, the PNs have been used to solve the voltage optimization problem in its reduced Equation (24). To this aim, 6 DERs are connected to the grid at nodes #6, #11, #13, #18, #21, #24, as shown in Figure 1b. Two DER configurations are considered to evaluate the impact of DER active and reactive powers injections/absorptions on the optimization of the voltage profiles: • DG-DERs are DGs connected to the grid by inverters; they inject a fixed active power =20 kW and are equipped with a Volt/Var control with =±15 kVAr (the European Standards [31] impose that inverters of significant powers (greater than 11 kW) have a "rectangular capability chart", which guarantees independence of the two powers in their admissible ranges. In the considered case, the inverter guarantees a "rectangular" capability chart that allows to vary the reactive power in the range [−15, 15] kVAr, independently from the value of the active power); in the initial operating point =0 kVAr. • DG and BESS-DERs are DGs connected to the grid by inverters that inject a fixed active power =10 kW and that are equipped with a Volt/Var control with = ±15 kVAr (in the initial operating point = 0 kVAr); DGs connected to nodes #6, #11, #18, #24 are also equipped with BESSs injecting/absorbing a controlled active power ; it is assumed that the state-ofcharge of the BESSs allows to vary in the range [−5, 5] kW.
To measure the performance of the Proposed VOP, the optimal solution must firstly be obtained. To this aim, for each Case and each DER configuration, firstly the minimization problem (23), optimizing all the nodal voltages with the PF equations as constraints, is solved by optimal power flow algorithm in MATPOWER [32]. In these problems, as well as in all the problems in the remainder, it is assumed , = 1.0 p.u. for all i. The optimal solution is expressed in terms of setpoints of the powers injected by DERs and the resulting voltage profile is the optimal one for the considered case, thus it is referred to as Benchmark. Then, the proposed problem (24), optimizing the voltages of the only PNs with the constraints derived from the linear model of the distribution system, is solved by classical quadratic programming algorithm in MATLAB ® ; in the following this problem is referred to as proposed voltage optimization problem (VOP).
Numerous VOPs have been formulated and solved considering different numbers of VCZs (from three to six) and different clustering algorithms, in particular the P and the Q algorithms as described in the previous Section 3.1. The solutions of the proposed VOP are the variations of the setpoints of the powers injected by DERs with respect to the initial operating point. Once the set-points are derived, such power injections are used to solve a power flow in MATPOWER so as to obtain the resulting voltage profile, which is also used to evaluate the value of the objective function of problem

Voltage Optimization Problem
To test the effectiveness of dividing the distribution network into VCZs, the PNs have been used to solve the voltage optimization problem in its reduced Equation (24). To this aim, 6 DERs are connected to the grid at nodes #6, #11, #13, #18, #21, #24, as shown in Figure 1b. Two DER configurations are considered to evaluate the impact of DER active and reactive powers injections/absorptions on the optimization of the voltage profiles: • DG-DERs are DGs connected to the grid by inverters; they inject a fixed active power P dg i = 20 kW and are equipped with a Volt/Var control with Q dg i = ±15 kVAr (the European Standards [31] impose that inverters of significant powers (greater than 11 kW) have a "rectangular capability chart", which guarantees independence of the two powers in their admissible ranges. In the considered case, the inverter guarantees a "rectangular" capability chart that allows to vary the reactive power in the range [−15, 15]  To measure the performance of the Proposed VOP, the optimal solution must firstly be obtained. To this aim, for each Case and each DER configuration, firstly the minimization problem (23), optimizing all the nodal voltages with the PF equations as constraints, is solved by optimal power flow algorithm in MATPOWER [32]. In these problems, as well as in all the problems in the remainder, it is assumed V re f ,i = 1.0 p.u. for all i. The optimal solution is expressed in terms of set-points of the powers injected by DERs and the resulting voltage profile is the optimal one for the considered case, thus it is referred to as Benchmark. Then, the proposed problem (24), optimizing the voltages of the only PNs with the constraints derived from the linear model of the distribution system, is solved by classical quadratic programming algorithm in MATLAB ® ; in the following this problem is referred to as proposed voltage optimization problem (VOP).
Numerous VOPs have been formulated and solved considering different numbers of VCZs (from three to six) and different clustering algorithms, in particular the P and the Q algorithms as described in the previous Section 3.1. The solutions of the proposed VOP are the variations of the set-points of the powers injected by DERs with respect to the initial operating point. Once the set-points are derived, such power injections are used to solve a power flow in MATPOWER so as to obtain the resulting voltage profile, which is also used to evaluate the value of the objective function of problem (23) in such suboptimal solution. Some results are presented in the following to evidence the main results that have been obtained.
Always referring to the Case R-DG, the nodal voltage profiles are shown, along feeder 1 for P and Q algorithms, respectively, in Figures 4 and 5, along feeder 2 for P and Q algorithms, respectively, in Figures 6 and 7. Different profiles are plotted in each graph: the blue line refers to the voltage profile with DGs and Q dg re f ,i = 0 for all i; the red line refers to the Benchmark; the other lines refer to the voltage profiles obtained by the proposed VOP applied to different numbers of VCZs (from 3 to 6). The blue profile is not strictly decreasing along both the feeders because of the presence of DGs. In particular, DGs connected to the feeder 1 determine an increase of the voltage amplitudes in the main from nodes #3 to #7 as well as in the lateral 1 at nodes #10 and #11 and in the lateral 2 at node #13. Similarly, DGs connected to the feeder 2 determine an increase of the voltage amplitudes from nodes #15 to #18 in the main, at nodes #20 and #21 in the lateral 1 and at nodes #23 and #24 in the lateral 2. Despite the voltage rise, the voltage amplitudes are within the acceptable range for all the nodes.
Energies 2018, 11, x FOR PEER REVIEW 16 of 28 (23) in such suboptimal solution. Some results are presented in the following to evidence the main results that have been obtained. Always referring to the Case R-DG, the nodal voltage profiles are shown, along feeder 1 for P and Q algorithms, respectively, in Figures 4 and 5, along feeder 2 for P and Q algorithms, respectively, in Figures 6 and 7. Different profiles are plotted in each graph: the blue line refers to the voltage profile with DGs and , = 0 for all i; the red line refers to the Benchmark; the other lines refer to the voltage profiles obtained by the proposed VOP applied to different numbers of VCZs (from 3 to 6). The blue profile is not strictly decreasing along both the feeders because of the presence of DGs. In particular, DGs connected to the feeder 1 determine an increase of the voltage amplitudes in the main from nodes #3 to #7 as well as in the lateral 1 at nodes #10 and #11 and in the lateral 2 at node #13. Similarly, DGs connected to the feeder 2 determine an increase of the voltage amplitudes from nodes #15 to #18 in the main, at nodes #20 and #21 in the lateral 1 and at nodes #23 and #24 in the lateral 2. Despite the voltage rise, the voltage amplitudes are within the acceptable range for all the nodes.     (23) in such suboptimal solution. Some results are presented in the following to evidence the main results that have been obtained. Always referring to the Case R-DG, the nodal voltage profiles are shown, along feeder 1 for P and Q algorithms, respectively, in Figures 4 and 5, along feeder 2 for P and Q algorithms, respectively, in Figures 6 and 7. Different profiles are plotted in each graph: the blue line refers to the voltage profile with DGs and , = 0 for all i; the red line refers to the Benchmark; the other lines refer to the voltage profiles obtained by the proposed VOP applied to different numbers of VCZs (from 3 to 6). The blue profile is not strictly decreasing along both the feeders because of the presence of DGs. In particular, DGs connected to the feeder 1 determine an increase of the voltage amplitudes in the main from nodes #3 to #7 as well as in the lateral 1 at nodes #10 and #11 and in the lateral 2 at node #13. Similarly, DGs connected to the feeder 2 determine an increase of the voltage amplitudes from nodes #15 to #18 in the main, at nodes #20 and #21 in the lateral 1 and at nodes #23 and #24 in the lateral 2. Despite the voltage rise, the voltage amplitudes are within the acceptable range for all the nodes.     Referring to the first DER configuration (DG) when the resistive parameters of the LV network are prevailing (Case R-DG), Table 8 reports the values of the objective function obtained by the proposed VOP and normalized with respect the Benchmark, which in this case is equal to 7.3 10 -3 p.u.. As evident from Table 8, the obtained values are of the same order of magnitude as the Benchmark: thus, the results can be considered substantially equivalent, with a slightly-better performance of the Q algorithm. It is interesting to underline that, in spite of the similarity of the values of the objective function, the proposed VOP provides solutions which significantly differ from the Benchmark in terms of the set points for the DG reactive powers , at nodes #11, #18, #21, as shown in Table 9.    Referring to the first DER configuration (DG) when the resistive parameters of the LV network are prevailing (Case R-DG), Table 8 reports the values of the objective function obtained by the proposed VOP and normalized with respect the Benchmark, which in this case is equal to 7.3 10 -3 p.u.. As evident from Table 8, the obtained values are of the same order of magnitude as the Benchmark: thus, the results can be considered substantially equivalent, with a slightly-better performance of the Q algorithm. It is interesting to underline that, in spite of the similarity of the values of the objective function, the proposed VOP provides solutions which significantly differ from the Benchmark in terms of the set points for the DG reactive powers , at nodes #11, #18, #21, as shown in Table 9.  Referring to the first DER configuration (DG) when the resistive parameters of the LV network are prevailing (Case R-DG), Table 8 reports the values of the objective function obtained by the proposed VOP and normalized with respect the Benchmark, which in this case is equal to 7.3 10 −3 p.u. As evident from Table 8, the obtained values are of the same order of magnitude as the Benchmark: thus, the results can be considered substantially equivalent, with a slightly-better performance of the Q algorithm. It is interesting to underline that, in spite of the similarity of the values of the objective function, the proposed VOP provides solutions which significantly differ from the Benchmark in terms of the set points for the DG reactive powers Q dg re f ,i at nodes #11, #18, #21, as shown in Table 9.   Concerning to the Case R-DG and BESS, Table 10 reports the obtained values of the objective function normalized to the Benchmark, which in this case is equal to 0.12·10 −3 p.u. and, then, significantly lower than the one in the previous Case R-DG. This is due to the possibility of varying also the active powers of BESSs which has a significant impact in a predominately-resistive network. From Table 10 it is evident that the results are substantially equivalent to the Benchmark with the only exceptions of the clustering into three and four VCZs provided by the P alghoritm which are significantly worse; for five VCZs the P algorithm shows a better result than the Q algorithm. Table 11 reports the resulting set-points for the active power aborbed/injected by BESS P bess re f ,i and for the reactive power absorbed/injected by DG Q dg re f ,i . The active powers injected at nodes #13 and #21 are always equal to zero since no BESSs are installed. The optimal solution of the Benchmark is different from the ones given by the proposed VOP; indeed, the set points P bess re f ,i are rather similar whereas the set points Q dg re f ,i of the proposed VOP hardly ever reach the saturation values differently from the Benchmark, as highlighted in bold in Table 11. Always referring to the Case R-DG and BESS, the nodal voltage profiles are shown, along feeder 1 for P and Q algorithms, respectively, in Figures 8 and 9, along feeder 2 for P and Q algorithms, respectively, in Figures 10 and 11. In absence of optimization (blue line) the voltage profile is within the acceptable range, although characterized by a voltage rise at the nodes with DGs. All the optimizations flatten the voltage profiles around the reference values equal to 1.0 p.u.             The optimal voltage profile (red line) is obtained by injecting reactive powers by DGs and by absorbing active powers by BESSs. In fact, the reactive power mainly affects the voltage drop of the MV/LV transformer; then, the reactive power injections allow to increase the voltage at the LV busbar of the substation around 1.0 p.u. As counter effect, the voltages would rise above 1.0 p.u. along the feeders where DGs are connected; then, the active power absorptions by BESSs flatten the nodal voltages around the reference values, since in a resistive network the impact of the BESS on the voltage drops along the feeders is dominant. The values of the nodal voltage amplitudes obtained by the proposed VOP are quite similar to the optimal ones, except for the cases of three VCZs and four VCZs clustered according to the P algorithm, in which the voltage profile of the feeder 2 is far from the reference values, thus explaining the results in Table 10. Similarly to the previous case, the voltage profiles improve when the number of the VCZs increases. Table 12 shows the objective functions of the Benchmark and the proposed VOP for, the Case R≈X with reference to the two DER configurations. The values are normalized with respect to the corresponding Benchmarks, which are equal to 2.1·10 −3 p.u. and 0.12·10 −3 p.u., respectively, for the configuration DG and the configuration DG and BESS. In a distribution network with comparable  The optimal voltage profile (red line) is obtained by injecting reactive powers by DGs and by absorbing active powers by BESSs. In fact, the reactive power mainly affects the voltage drop of the MV/LV transformer; then, the reactive power injections allow to increase the voltage at the LV busbar of the substation around 1.0 p.u. As counter effect, the voltages would rise above 1.0 p.u. along the feeders where DGs are connected; then, the active power absorptions by BESSs flatten the nodal voltages around the reference values, since in a resistive network the impact of the BESS on the voltage drops along the feeders is dominant. The values of the nodal voltage amplitudes obtained by the proposed VOP are quite similar to the optimal ones, except for the cases of three VCZs and four VCZs clustered according to the P algorithm, in which the voltage profile of the feeder 2 is far from the reference values, thus explaining the results in Table 10. Similarly to the previous case, the voltage profiles improve when the number of the VCZs increases. Table 12 shows the objective functions of the Benchmark and the proposed VOP for, the Case R≈X with reference to the two DER configurations. The values are normalized with respect to the corresponding Benchmarks, which are equal to 2.1·10 −3 p.u. and 0.12·10 −3 p.u., respectively, for the configuration DG and the configuration DG and BESS. In a distribution network with comparable The optimal voltage profile (red line) is obtained by injecting reactive powers by DGs and by absorbing active powers by BESSs. In fact, the reactive power mainly affects the voltage drop of the MV/LV transformer; then, the reactive power injections allow to increase the voltage at the LV busbar of the substation around 1.0 p.u. As counter effect, the voltages would rise above 1.0 p.u. along the feeders where DGs are connected; then, the active power absorptions by BESSs flatten the nodal voltages around the reference values, since in a resistive network the impact of the BESS on the voltage drops along the feeders is dominant. The values of the nodal voltage amplitudes obtained by the proposed VOP are quite similar to the optimal ones, except for the cases of three VCZs and four VCZs clustered according to the P algorithm, in which the voltage profile of the feeder 2 is far from the reference values, thus explaining the results in Table 10. Similarly to the previous case, the voltage profiles improve when the number of the VCZs increases. Table 12 shows the objective functions of the Benchmark and the proposed VOP for, the Case R≈X with reference to the two DER configurations. The values are normalized with respect to the corresponding Benchmarks, which are equal to 2.1·10 −3 p.u. and 0.12·10 −3 p.u., respectively, for the configuration DG and the configuration DG and BESS. In a distribution network with comparable resistive and inductive line parameters, the evaluation of VCZs resulting from the application of the P and Q algorithms are equivalent in terms of the resulting values for the objective function. In fact, for both the DER configurations, the two zoning algorithms yield the same solution, except for the case of four VCZs. However, the errors of the proposed VOP become acceptable starting from five VCZs. As an example, Figures 12 and 13 show the nodal voltage profiles of, respectively, feeder 1 and feeder 2 for the benchmark and the proposed VOP with four VCZs in the Case R≈X-DG. From Figure 12 the overlapping of the voltage profiles along feeder 1 obtained by the proposed VOP for the P and the Q algorithms is evident. Conversely, along the feeder 2 the benchmark is well followed by the P algorithm as far as the main and the lateral 1 are considered, whereas the Q algorithm presents a better performance for the lateral 2. The poor performance of the proposed VOP with the P algorithm is related to the absence of a PN belonging to lateral 2 (see Table 5) and is confirmed by the high value of the objective function (see Table 12).  Table 13 shows the normalized objective functions of the Benchmark and the proposed VOP for the Case X with reference to the two DER configurations. In general, the results obtained using the Q algorithm are better than the ones of the P algorithm, as expectable in a predominantly-inductive distribution network. The values of the objective function for the P algorithm are of one order of magnitude greater than the Benchmark, except for the case of three VCZs and for the case of six VCZs. Figures 14 and 15 show the nodal voltage profiles of, respectively, feeder 1 and feeder 2 for the Benchmark and the Proposed VOP with three VCZs in the Case X-DG. The P algorithm provides the best approximation of the optimal solution when three VCZs are adopted and this result is coherent with the better performance of Silhouette index for the P algorithm for 3 VCZs in the Case X, see Figure 3. This is due to the choice of node #20 as PN by the Q algorithm that gives rise to a high difference of the voltage amplitudes at node #21 of the Benchmark and the proposed VOP ( Figure 14). From Table 13, it would seem that results have become acceptable starting from six VCZs. Actually, taking into account that the objective functions of the Benchmark are equal to 0.50 10 −3 p.u. and 0.08 10 −3 p.u. for, respectively, the DGs and the DGs with BESS configurations, the voltage regulation provides good performance also for a reduced number of VCZs. As an example, for the Case X-DG and BESS the nodal voltage profiles are shown, along feeder 1 for P and Q algorithms, respectively, in Figures 16  and 17, along feeder 2 for P and Q algorithms, respectively, in Figures 18 and 19.

Conclusions
The problem of the voltage profile optimization in a distribution system including distributed energy resources has been tackled adopting a centralized approach. To avoid the monitoring of a large number of electrical quantities in the distribution network as well as the solution of a large nonlinear programming problem, the voltage optimization has been re-formulated as a quadratic programming of reduced dimension based on control zones. The voltage control zones have been identified on the basis of the structural characteristics of the power grid (independently from the number and position of the distributed energy resources) and account for the sensitivity of the nodal voltages on active and reactive powers. General considerations can be made, although derived from a case study. First of all, the proposed voltage optimization problem always guarantees an improvement of the voltage profile; it gives a suboptimal solution which is quite similar to the optimal one provided that an adequate number of voltage control zones is used. Concerning the choice of the algorithm for the aggregation of the voltage control zones, the use of the voltage sensitivities on the reactive power is preferable if the voltage control acts on the reactive power injections/absorptions, as easily foreseeable. When the control acts on both active and reactive powers, the choice is to be related to the R/X ratio of the distribution lines. In fact, in predominantly-

Conclusions
The problem of the voltage profile optimization in a distribution system including distributed energy resources has been tackled adopting a centralized approach. To avoid the monitoring of a large number of electrical quantities in the distribution network as well as the solution of a large nonlinear programming problem, the voltage optimization has been re-formulated as a quadratic programming of reduced dimension based on control zones. The voltage control zones have been identified on the basis of the structural characteristics of the power grid (independently from the number and position of the distributed energy resources) and account for the sensitivity of the nodal voltages on active and reactive powers. General considerations can be made, although derived from a case study. First of all, the proposed voltage optimization problem always guarantees an improvement of the voltage profile; it gives a suboptimal solution which is quite similar to the optimal one provided that an adequate number of voltage control zones is used. Concerning the choice of the algorithm for the aggregation of the voltage control zones, the use of the voltage sensitivities on the reactive power is preferable if the voltage control acts on the reactive power injections/absorptions, as easily foreseeable. When the control acts on both active and reactive powers, the choice is to be related to the R/X ratio of the distribution lines. In fact, in predominantly-

Conclusions
The problem of the voltage profile optimization in a distribution system including distributed energy resources has been tackled adopting a centralized approach. To avoid the monitoring of a large number of electrical quantities in the distribution network as well as the solution of a large nonlinear programming problem, the voltage optimization has been re-formulated as a quadratic programming of reduced dimension based on control zones. The voltage control zones have been identified on the basis of the structural characteristics of the power grid (independently from the number and position of the distributed energy resources) and account for the sensitivity of the nodal voltages on active and reactive powers. General considerations can be made, although derived from a case study. First of all, the proposed voltage optimization problem always guarantees an improvement of the voltage profile; it gives a suboptimal solution which is quite similar to the optimal one provided that an adequate number of voltage control zones is used. Concerning the choice of the algorithm for the aggregation of the voltage control zones, the use of the voltage sensitivities on the reactive power is preferable if the voltage control acts on the reactive power injections/absorptions, as easily foreseeable. When the control acts on both active and reactive powers, the choice is to be related to the R/X ratio of the distribution lines. In fact, in predominantly-inductive networks, such as MV networks, better performance is achieved if the aggregation is based on sensitivities to reactive powers. On the other hand, in the case of LV distribution networks with predominantly-resistive lines, the distributed Energies 2019, 12, 390 26 of 28 energy resources impact on the substation busbar voltage mainly through the reactive powers and on the shape of the voltage profile along the feeders mainly through the active powers; consequently, it is advisable to use an aggregation algorithm based on the sensitivities to both active and reactive powers. Finally, the proposed formulation of the voltage optimization problem based on voltage control zones is by its nature oriented to the decentralized control, whose implementation is the future step of this research.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The main steps of the extended hierarchical algorithm are: i.
Input: D P and D Q ; ii.
Initialize j = 1 and D P (j) = D P and D Q (j) = D Q ; iii.
Find the minimum electrical distance D min m,n = 0 in both D P (j) and D Q (j); iv.
Aggregate in D P (j) and D Q (j) the nodes m and n in a fictious node by substituting D P m,n (j) = 0 and D Q m,n (j) = 0 and by imposing symmetrical distances D P n,m (j) = D P m,n (j) and D Q n,m (j) = D Q m,n (j); v.
Re-compute firstly in D P (j) and then in D Q (j) the electrical distance different from zero among the fictious nodes m and n and the remaining nodes i = 1, . . . , N with the maximum distance criterion (the adoption of the maximum distance criterion assures that, in the step iii. of the iteration j, the minimum distance is larger than the one evaluated in the step iii. of the iteration j − 1) and impose symmetrical distances. As an example for D P (j): FOR i = 1:N IF D P i,m (j) = 0 AND D P i,n (j) = 0 THEN D P i,m (j) = max D P i,m (j), D P i,n (j) D P i,n (j) = D P i,m (j) END D P m,i (j) = D P i,m (j) D P n,i (j) = D P i,n (j) END vi. Store D P (j) and D Q (j), representative of a network divided in N − 1 VCZs; vii.
If D P (j) and D Q (j) are not null matrices go back to step iii. with j = j + 1; viii.
Output: N matrices D P (j) and D Q (j) representative of a network divided in N-1, . . . ,1 VCZs; ix.
Fix the number N vcz of VCZs on the basis of clustering quality indices or of the type of control problem; x.
If the i-th row corresponds to a node which is already included in a VCZ PQ k for k < h, then go to step xi. xiii.
Insert in the h-th control area VCZ PQ h the i-th node and all the nodes k corresponding to the columns of D P (N vcz ) such that D P i,k (j) = 0 for k = 1, . . . , N (the same result is obtained by using D Q (N VCZ ) in place of D P (N VCZ ));