Sensitivity-Based Model of Low Voltage Distribution Systems with Distributed Energy Resources

Abstract

A key issue in Low Voltage (LV) distribution systems is to identify strategies for the optimal management and control in the presence of Distributed Energy Resources (DERs). To reduce the number of variables to be monitored and controlled, virtual levels of aggregation, called Virtual Microgrids (VMs), are introduced and identified by using new models of the distribution system. To this aim, this paper, revisiting and improving the approach outlined in a conference paper, presents a sensitivity-based model of an LV distribution system, supplied by an Medium/Low Voltage (MV/LV) substation and composed by several feeders, which is suitable for the optimal management and control of the grid and for VM definition. The main features of the proposed method are: it evaluates the sensitivity coefficients in a closed form; it provides an overview of the sensitivity of the network to the variations of each DER connected to the grid; and it presents a limited computational burden. A comparison of the proposed method with both the exact load flow solutions and a perturb-and-observe method is discussed in a case study. Finally, the method is used to evaluate the impact of the DERs on the nodal voltages of the network.

1. Introduction

Existing distribution networks present an inflexible structure which makes them inadequate to allow a large spread of Distributed Energy Resources (DERs), such as distributed generators from renewable energy sources, cogeneration facilities, energy storage systems (including electric vehicles), and controllable loads (i.e., adopting active demand) [1,2,3]. To overcome this problem, the smart grid paradigm is asserted and promoted by extensive use of Information and Communication Technologies [4,5,6]. However, the advent of smart distribution grids will not take place by a revolution, because it would require expensive investments, especially for Low Voltage (LV) networks that actually present a low level of control and automation. The evolution from the present networks to the future smart grids will be progressive instead. Then, the main technical challenge is to improve the flexibility of the existing LV networks with limited investments so as to allow the wide spread of DERs and enhance the power quality of the grids [7].

In this context, a key issue is to identify strategies for the optimal management and control of the LV distribution networks in the presence of DERs. For an optimal distribution system operation, the same strategy adopted at the transmission level could be used, in which a single operator controls the whole power system by gathering measurements and acting on each electrical node of the grid. Unfortunately, if such an approach were applied to LV distribution networks, it would be technically and economically infeasible because of the need to acquire and process a huge amount of information, due to the large number of nodes and customers in the grids.

To reduce the number of variables to be monitored and controlled, an alternative approach consists of modeling each LV distribution system by a cluster of Virtual Microgrids (VMs) and building communication infrastructures among a few distributed controllers of the VMs and a control center located at the Medium/Low Voltage (MV/LV) substation. Similarly to a microgrid [8], a VM is composed of parts of the distribution networks including different types of DERs and their control systems. However, differently from a microgrid, a VM is not self sustained from an energy point of view, and it is not equipped with an autonomous centralized management and control system. To identify VMs, new models of the distribution system are required. These models must be able to overcome the classical formulation adopted at the transmission level (based on non-linear equations written for each node of the grid) so as to facilitate the aggregation of portions of network on the basis of the impact of DERs on specific control actions (i.e., voltage regulation [9,10] and current flow limitation [11]). Sensitivity analysis can help face such a problem by quantifying the dependency between nodal voltages and power flows of the grid as a linear function of both load consumptions and DERs’ absorbtions/injections.

Several methods have been proposed in literature to perform sensitivity analysis. In [9], a Jacobian-based method has been shown to control the nodal voltages by regulating the reactive powers injected by distributed generators. Voltage magnitude as well as active and reactive power-loss sensitivity coefficients for radial distribution networks have been calculated in [11] by using only a base-case load flow solution. In [12], both the Gauss–Seidel method and Z-bus matrix have been used to derive voltage and loss sensitivity coefficients as functions of real and reactive power injections. In [13], voltage sensitivity factors and a linear load flow model for LV distribution system have been derived by using historical smart meter data. In [14], sensitivities of bus voltage magnitudes to load changes have been calculated using the adjoint-network method. In [15], a perturb-and-observe algorithm has been proposed to compensate voltage rises induced by distributed generators and to test siting and sizing of capacitors. Initial studies on the application of a closed-form solution of the Dist–Flow equations of an LV distribution system have been shown in [16,17].

This paper, adopting the approach outlined in [16], proposes a sensitivity-based model of an LV distribution system supplied by an MV/LV substation and composed by several feeders, which is suitable for the optimal management and control of the grid and, then, for VM definition. Such a model provides, in a closed form, the variations of the electrical variables characterizing each node of a feeder (e.g., out-flowing active and reactive powers, nodal voltage) as linear functions of the powers injected/absorbed by each DER connected to the LV network. The key steps of the proposed model are: (i) the linearization of the branch flow equations of the electrical components of the distribution system; (ii) the use of the chain-rule; (iii) the use of the border constraints at both the MV substation busbar and the end-nodes of the feeders; and (iv) the use of the coupling equations at the LV substation busbar. The main features of the proposed method are the following: (i) it evaluates the sensitivity coefficients in a closed form; (ii) it provides an overview of the sensitivity of the network to the variations of each DER connected to the grid; and (iii) it presents a limited computational burden. The last feature is very important because of the large number of nodes in LV distribution networks, of the need to repeat the analysis when the operating conditions of the distribution system change, and of the large dimension of the optimal management problems that account for the random nature of DERs.

With respect to [16], the paper presents in more details the method derived from the model, puts it into the frame of the sensitivity analysis, compares it with other methods, and applies it to DER impact evalutation in view of VM definition. After a brief introduction, in the paper, the sensitivity-based model of the LV distribution system is recalled and the step-by-step procedure of the method is presented in Section 2; then, an extensive case study is discussed in Section 3. Firstly, different sensitivity matrices are evaluated to illustrate how the proposed model provides the response of the LV distribution system to the DERs’ injections. Then, a comparison of the proposed method with both exact load flow solution and a perturb-and-observe method is discussed. Finally, the proposed method is used to evaluate the impact of the DERs on the nodal voltages of the network so as to provide a road map to perform voltage regulation. Conclusions and further research developments are reported in Section 4.

2. Distribution System Modeling

Let us consider the LV distribution system reported in Figure 1: an MV/LV substation supplies m feeders; the h-th feeder is composed of $n_h$ branches and $n_h+1$ nodes (node 0 is common to all the feeders being the LV busbar of the substation); along the feeders, uncontrolled loads and/or DERs can be connected. In balanced operating conditions, the MV/LV supplying system is modeled by the electric circuit shown in Figure 2. Two electric parameters identify the supplying system, namely the resistance $R_{tr}$ of the MV/LV transformer and the equivalent impedance $X_{eq}=X_{tr}+X_{sc}$, where $X_{tr}$ is the transformer reactance and $X_{sc}$ the short-circuit impedance of the MV node. Three electric variables characterize the model at the supplying MV node, namely the open-circuit voltage amplitude $V_{MV}$ of the MV node (assumed to be fixed as the slack node), the in-flowing active power $P_{MV}$ and the in-flowing reactive power $Q_{MV}$. Analogously, three electric variables characterize the MV/LV substation at the receiving node 0, namely the voltage amplitude $V_{tr}$ and the out-flowing active and reactive powers, $P_{tr}$ and $Q_{tr}$, respectively. In a similar way, the j-th branch of a generic feeder of the LV distribution system is represented by the electric circuit shown in Figure 3. Two electric parameters identify the branch, namely the resistance $R_j$ and the reactance $X_j$; three electric variables characterize the branch at the supplying node $j-1$, namely the voltage amplitude $V_{j-1}$, the in-flowing active power $P_{j-1}$ and the in-flowing reactive power $Q_{j-1}$; analogously, the same electric variables characterize the branch at the receiving node j, namely the voltage $V_j$, and the out-flowing active and reactive powers, $P_j$ and $Q_j$, respectively; in addition, other two external variables are taken into account at the receiving node j, which are the shunt active power $P_{S,j}$ and the shunt reactive power $Q_{S,j}$ defined as:

$$\begin{gathered} P_{S,j}=P_{DER,j}-P_{L,j}, \\ {Q}_{S,j}=Q_{DER,j}-Q_{L,j}, \end{gathered}$$

where $P_{DER,j}$ and $Q_{DER,j}$ are, respectively, the active and reactive powers injected by the DER connected to the j-th node; $P_{L,j}$ and $Q_{L,j}$ are, respectively, the active and reactive powers absorbed by the uncontrolled load connected to the j-th node, which are assumed to be constant. The model of the LV distribution system, which is useful for the VM definition and which overcomes the limits of the classical formulation adopted at the transmission level (based on non-linear equations written for each node of the grid), is presented in the following. In particular, by the linearization of the branch flow equations associated with the components of the LV distribution system represented in Figure 2 and Figure 3, and by imposing both border and coupling conditions, a closed form solution of the linearized DistFlow equations is obtained through the use of the chain rule. Such a closed form represents the linear model of the distribution system with DERs, and it is expressed in terms of the variations of the electrical variables of the LV networks (with respect to an initial operating point) as linear functions of the DER injections. The model is built step by step: firstly, the model of a generic branch and of a feeder composed of different branches is developed; then, the model of the MV/LV supplying system is introduced and combined with the models of the feeders, yielding the model of the whole LV distribution system.

2.1. Branch Model

The model of the j-th branch is derived from the following steady-state branch flow equations [18]:

$$\begin{gathered} P_j=P_{j-1}-R_j(P_{j-1}^2+Q_{j-1}^2)/V_{j-1}^2+P_{S,j}, \\ {Q}_j=Q_{j-1}-X_j(P_{j-1}^2+Q_{j-1}^2)/V_{j-1}^2+Q_{S,j}, \\ \begin{array}{ll}{V}_j^2& =V_{j-1}^2-2\left(R_j{P}_{j-1}+X_j{Q}_{j-1}\right)\\ & +(R_j^2+X_j^2)(P_{j-1}^2+Q_{j-1}^2)/V_{j-1}^2.\end{array} \end{gathered}$$

(1)

Since they are not linear, they can be linearized around an initial operating point of the distribution system, characterized by null powers injected by all the DERs connected to the LV distribution system and a given absorption of all the uncontrolled loads. In the following, the electrical variables evaluated in the initial operating point are identified by zero superscript.

Let the vector of nodal variables ${\mathit{x}}_j$ be defined as:

$${\mathit{x}}_j={\left(P_j{Q}_j{V}_j^2\right)}^{\mathrm{T}}.$$

The linearization of Equation (1) yields the following branch model:

$$\Delta {\mathit{x}}_j={\mathbf{J}}_j\Delta {\mathit{x}}_{j-1}+{(\Delta P_{S,j}\Delta Q_{S,j}0)}^{\mathrm{T}},$$

(2)

$\Delta {\mathit{x}}_j$ being the variation of the nodal variables ${\mathit{x}}_j$ with respect to their value ${\mathit{x}}_j^0$ in the initial operating point

$$\Delta {\mathit{x}}_j={\mathit{x}}_j-{\mathit{x}}_j^0={\left(\Delta P_j\Delta Q_j\Delta V_j^2\right)}^{\mathrm{T}},$$

with ${\mathbf{J}}_j$ being the Jacobian matrix related to (1) evaluated in the initial operating point and $\Delta P_{S,j}$, $\Delta Q_{S,j}$ the injections of the active and reactive powers by the DERs connected at the receiving node.

Equation (2) is the branch model: it expresses the vector $\Delta {\mathit{x}}_j$ of the variation of the electrical variables at the receiving node ${\mathit{x}}_j$ with respect to the initial operating point ${\mathit{x}}_j^0$ as a function of

• the vector $\Delta {\mathit{x}}_{j-1}$ of the variation of the electrical variables at the supplying node ${\mathit{x}}_{j-1}$, with respect to the initial operating point ${\mathit{x}}_{j-1}^0$, through the Jacobian matrix ${\mathbf{J}}_j$; and
• the vector ${(\Delta P_{S,j}\Delta Q_{S,j}0)}^{\mathrm{T}}$ of the DERs active and reactive power injections at the receiving node, which are assumed to be assigned enforcements.

2.2. Feeder Model

The electric equivalent circuit of the h-th feeder of the LV distribution system is composed of a series of $n_h$ branches (Figure 4). Then, the model of the feeder is composed of $n_h$ equations of type (2):

$$\Delta {\mathit{x}}_j={\mathbf{J}}_j\Delta {\mathit{x}}_{j-1}+{\left(\Delta P_{S,j}\Delta Q_{S,j}0\right)}^{\mathrm{T}},$$

(3)

with $j=1,...,n_h$.

Applying the chain rule to Equation (3),

$$\Delta {\mathit{x}}_j={\mathbf{M}}_j\Delta {\mathit{x}}_0+\sum _{k=1}^j{\mathbf{N}}_{j,k}{(\Delta P_{S,k}\Delta Q_{S,k}0)}^{\mathrm{T}}$$

(4)

is induced where

$${\mathbf{M}}_j=\prod _{\ell =0}^{j-1}{\mathbf{J}}_{j-\ell }j=1,\dots ,n_h,$$

(5)

$${\mathbf{N}}_{j,k}=\prod _{\ell =0}^{j-k-1}{\mathbf{J}}_{j-\ell }j=1,\dots ,n_{h}k=1,\dots ,j-1,$$

(6)

and ${\mathbf{N}}_{j,j}$ is equal to the identity matrix I. Assuming ${\mathbf{N}}_{j,k}=0$ for $k>j$, Equation (4) can be re-written as

$$\Delta {\mathit{x}}_j={\mathbf{M}}_j\Delta {\mathit{x}}_0+\sum _{k=1}^{n_h}{\mathbf{N}}_{j,k}\left(\begin{array}{c}\Delta P_{S,k}\\ \Delta Q_{S,k}\\ 0\end{array}\right),$$

(7)

with $j=1,\dots ,n_h$.

Equation (7) is a set of $n_h$ equations in $n_h+1$ variables, which are $\Delta {\mathit{x}}_j$ and $\Delta {\mathit{x}}_0$. To obtain a defined problem, three additional conditions must be assigned. They can be derived from the behavior of the feeder at its borders. In particular, the voltage amplitude at the head of the feeder is imposed by the MV/LV supplying substation and the active and reactive powers flowing out of the end-node of the feeder are always null. The border conditions can be expressed in terms of variations of the electric nodal variables as

$$\Delta V_0^2=\Delta {\tilde{V}}_0^2,\Delta P_{n_h}=0,\Delta Q_{n_h}=0,$$

(8)

where $\Delta {\tilde{V}}_0^2$ is the variation of the squared voltage amplitude imposed at the LV busbar by the MV/LV supplying substation. By substituting Equation (8) into Equation (7) for $j=n_h$,

$$\left(\begin{array}{c}0\\ 0\\ \Delta V_{n_h}^2\end{array}\right){=}_{{\mathbf{M}}_{n_h}}\left(\begin{array}{c}\Delta P_0\\ \Delta Q_0\\ \Delta {\tilde{V}}_0^2\end{array}\right)+\sum _{k=1}^{n_h}{\mathbf{N}}_{n_h,k}\left(\begin{array}{c}\Delta P_{S,k}\\ \Delta Q_{S,k}\\ 0\end{array}\right)$$

(9)

is obtained.

The subset composed of the first two rows in Equation (9) can be solved in the two unknowns, $\Delta P_0$ and $\Delta Q_0$, yielding:

$$\left(\begin{array}{c}\Delta P_0\\ \Delta Q_0\end{array}\right)=\sum _{k=1}^{n_h}{\overline{\mathbf{A}}}_k\left(\begin{array}{c}\Delta P_{S,k}\\ \Delta Q_{S,k}\end{array}\right)+\overline{\mathit{a}}\Delta {\tilde{V}}_0^2,$$

(10)

where,

$${\overline{\mathbf{A}}}_k=-{\overline{\mathbf{M}}}_{n_h}^{-1}{\overline{\mathbf{N}}}_{n_h,k},\overline{\mathit{a}}=-{\overline{\mathbf{M}}}_{n_h}^{-1}{\overline{\mathbf{m}}}_{n_h},$$

(11)

with ${\overline{\mathbf{M}}}_{n_h}$ and ${\overline{\mathbf{N}}}_{n_h,k}$ the $(2\times 2)$ leading principal minors of ${\mathbf{M}}_{n_h}$ and ${\mathbf{N}}_{n_h,k}$, respectively, and with ${\overline{\mathbf{m}}}_{n_h}$ the $(2\times 1)$ column vector composed of the first two rows in the third column of ${\mathbf{M}}_{n_h}$.

Combining Equation (10) with the condition on $\Delta V_0^2$ in Equation (8), the variable $\Delta {\mathit{x}}_0$ can be written in the form:

$$\Delta {\mathit{x}}_0=\left(\begin{array}{c}\Delta P_0\\ \Delta Q_0\\ \Delta V_0^2\end{array}\right)=\sum _{k=1}^{n_h}{\mathbf{A}}_k\left(\begin{array}{c}\Delta P_{S,k}\\ \Delta Q_{S,k}\\ 0\end{array}\right)+\mathit{a}\Delta {\tilde{V}}_0^2,$$

(12)

where each ${\mathbf{A}}_k$ is a $(3\times 3)$ matrix obtained by edging ${\overline{\mathbf{A}}}_k$ with a third row and column of zero elements, and $\mathit{a}$ is a $(3\times 1)$ vector obtained by adding a third row with a unitary element to $\overline{\mathit{a}}$.

By substituting $\Delta {\mathit{x}}_0$ given by Equation (12) into Equation (7) yields the model of the feeder:

$$\Delta {\mathit{x}}_j=\sum _{k=1}^{n_h}{\mathbf{B}}_{j,k}\left(\begin{array}{c}\Delta P_{S,k}\\ \Delta Q_{S,k}\\ 0\end{array}\right)+{\mathit{b}}_j\Delta {\tilde{V}}_0^2,$$

(13)

with $j=1,\dots ,n_h$, where:

$${\mathbf{B}}_{j,k}={\mathbf{M}}_j{\mathbf{A}}_k+{\mathbf{N}}_{j,k},{\mathit{b}}_j={\mathbf{M}}_j\mathit{a},$$

(14)

which are, respectively, ($3\times 3$) matrices and ($3\times 1$) vectors of sensitivity with known coefficients.

In Equation (13), the vectors $\Delta {\mathit{x}}_j$ of the variation of the electric variables at the receiving node of each branch belonging to the feeder are expressed as a function of:

• The vectors ${(\Delta P_{S,k}\Delta Q_{S,k}0)}^{\mathrm{T}}$ of the DERs active and reactive power injections in all the nodes of the feeder through the matrices ${\mathbf{B}}_{j,k}$; and
• The variation $\Delta {\tilde{V}}_0^2$ of the squared voltage amplitude at the LV busbar of the supplying substation through the vectors ${\mathit{b}}_j$.

To specify the model Equations (12) and (13) for each one of the m feeders, a superscript referring to the feeder number is added to the variables, matrices and vectors. Then, the model Equations (12) and (13) for the h-th feeders becomes:

$$\Delta {\mathit{x}}_0^h=\sum _{k=1}^{n_h}{\mathbf{A}}_k^h\left(\begin{array}{c}\Delta P_{S,k}^h\\ \Delta Q_{S,k}^h\\ 0\end{array}\right)+{\mathit{a}}^h\Delta {\tilde{V}}_0^2,$$

(15)

$$\Delta {\mathit{x}}_j^h=\sum _{k=1}^{n_h}{\mathbf{B}}_{j,k}^h\left(\begin{array}{c}\Delta P_{S,k}^h\\ \Delta Q_{S,k}^h\\ 0\end{array}\right)+{\mathit{b}}_j^h\Delta {\tilde{V}}_0^2,$$

(16)

with $j=1,\dots ,n_h$ and $h=1,...,m$.

2.3. Medium/Low Voltage Supplying System

Similarly to the branch model Equation (2), the MV/LV supplying system, can be modeled as:

$$\Delta {\mathit{x}}_{MV}={\mathbf{J}}^{sub}\Delta {\mathit{x}}_{tr},$$

(17)

where:

$$\begin{gathered} \Delta {\mathit{x}}_{MV}={(\Delta P_{MV}\Delta Q_{MV}0)}^{\mathrm{T}}, \\ \Delta {\mathit{x}}_{tr}={\left(\Delta P_{tr}\Delta Q_{tr}\Delta V_{tr}^2\right)}^{\mathrm{T}}, \end{gathered}$$

where $\Delta V_{MV}^2=0$ and ${\mathbf{J}}^{sub}$ is the Jacobian matrix of the linearized branch flow equations of the electric circuit in Figure 2.

2.4. Low Voltage Distribution System Model

To obtain the model of the whole LV distribution system in Figure 1, it is necessary to combine the model Equation of the MV/LV supplying system (17) with the model Equations of the feeders (15) and (16).

As evident from Figure 1, three coupling equations can be written:

$$\Delta V_{tr}^2=\Delta {\tilde{V}}_0^2;$$

(18)

$$\Delta P_{tr}=\sum _{i=1}^m\Delta P_0^i\Delta Q_{tr}=\sum _{i=1}^m\Delta Q_0^i.$$

(19)

By substituting Equations (18) and (19) into (17),

$$\left(\begin{array}{c}\Delta P_{MV}\\ \Delta Q_{MV}\\ 0\end{array}\right)={\mathbf{J}}^{sub}\left(\begin{array}{c}{\sum }_{i=1}^m\Delta P_0^i\\ {\sum }_{i=1}^m\Delta Q_0^i\\ \Delta {\tilde{V}}_0^2\end{array}\right)$$

(20)

is obtained.

Substituting the expression (10) generalized for m feeders, the third equation in (20) is written as:

$$0=J_{3,3}^{sub}\Delta {\tilde{V}}_0^2+\left(\begin{array}{c}{J}_{3,1}^{sub}{J}_{3,2}^{sub}\end{array}\right)\sum _{i=1}^m\left(\sum _{k=1}^{n_i}{\overline{\mathbf{A}}}_k^i\left(\begin{array}{c}\Delta P_{S,k}^i\\ \Delta Q_{S,k}^i\end{array}\right)+{\overline{\mathit{a}}}^i\Delta {\tilde{V}}_0^2\right),$$

(21)

with $(J_{3,1}^{sub},J_{3,2}^{sub},J_{3,3}^{sub})$ being the third row of ${\mathbf{J}}^{sub}$. By solving Equation (21) in the variable $\Delta {\tilde{V}}_0^2$,

$$\Delta {\tilde{V}}_0^2=-{\overline{\mathbf{\alpha }}}^{\mathbf{T}}\sum _{i=1}^m\sum _{k=1}^{n_i}{\overline{\mathbf{A}}}_k^i\left(\begin{array}{c}\Delta P_{S,k}^i\\ \Delta Q_{S,k}^i\end{array}\right)$$

(22)

is obtained, where $\overline{\mathbf{\alpha }}$ is a $(2\times 1)$ vector equal to:

$$\overline{\mathbf{\alpha }}={\left(\left(J_{3,1}^{sub}{J}_{3,2}^{sub}\right)\sum _{i=1}^m{\overline{\mathit{a}}}^i+J_{3,3}^{sub}\right)}^{-1}{\left(J_{3,1}^{sub}{J}_{3,2}^{sub}\right)}^{\mathrm{T}}.$$

(23)

By substituting Equation (22) into (16), the model of the whole LV distribution system is obtained:

$$\Delta {\mathit{x}}_j^h=\sum _{k=1}^{n_h}{\mathbf{B}}_{j,k}^h\left(\begin{array}{c}\Delta P_{S,k}^h\\ \Delta Q_{S,k}^h\\ 0\end{array}\right)+{\mathit{b}}_j^h\left(-{\overline{\mathbf{\alpha }}}^{\mathbf{T}}\sum _{i=1}^m\sum _{k=1}^{n_i}{\overline{\mathbf{A}}}_k^i\left(\begin{array}{c}\Delta P_{S,k}^i\\ \Delta Q_{S,k}^i\end{array}\right)\right),$$

(24)

with $j=1,\dots ,n_h$ and $h=1,...,m$.

By substituting into Equation (24) the matrix ${\mathbf{A}}_k^i$ from Equation (15) in place of ${\overline{\mathbf{A}}}_k^i$ and adding a third row with a null element to the vector ${(\Delta P_{S,k}^i\Delta Q_{S,k}^i)}^{\mathrm{T}}$ and to the vector $\overline{\mathbf{\alpha }}$ which becomes $\mathbf{\alpha }$, the model of the LV distribution system with DERs becomes:

$$\Delta {\mathit{x}}_j^h=\sum _{k=1}^{n_h}{\mathit{U}}_{j,k}^h\left(\begin{array}{c}\Delta P_{S,k}^h\\ \Delta Q_{S,k}^h\\ 0\end{array}\right)+\sum _{i=1,i\ne h}^m\sum _{k=1}^{n_i}{\mathit{T}}_{j,k}^{h,i}\left(\begin{array}{c}\Delta P_{S,k}^i\\ \Delta Q_{S,k}^i\\ 0\end{array}\right),$$

(25)

for $j=1,\dots ,n_h$ and $h=1,...,m$, where:

$${\mathit{U}}_{j,k}^h={\mathbf{B}}_{j,k}^h-{\mathit{b}}_j^h{\mathbf{\alpha }}^{\mathbf{T}}{\mathbf{A}}_k^h{\mathit{T}}_{j,k}^{h,i}=-{\mathit{b}}_j^h{\mathbf{\alpha }}^{\mathbf{T}}{\mathbf{A}}_k^i,$$

(26)

which are $(3\times 3)$ sensitivity matrices with known coefficients.

In Equation (25), the vectors $\Delta {\mathit{x}}_j^h$ of the variations of the electrical variables at the receiving node of the j-th branch in the h-th feeder are expressed as functions of:

• The vectors ${(\Delta P_{S,k}^h\Delta Q_{S,k}^{h}0)}^{\mathrm{T}}$ of the DERs active and reactive power injections in all the nodes of the h-th feeder through the matrices ${\mathit{U}}_{j,k}^h$; and
• The vectors ${(\Delta P_{S,k}^i\Delta Q_{S,k}^{i}0)}^{\mathrm{T}}$ of the DERs active and reactive power injections in all the nodes of the other $m-1$ feeders through the matrices ${\mathit{T}}_{j,k}^{h,i}$.

The second term in Equation (25) is due to the coupling among the feeders at the LV busbar of the MV/LV substation through the nodal voltage $V_{tr}$; a variation of the DER injections in any of the other $m-1$ feeders causes a variation of voltage at the supplying substation, and, in turn, of the electrical variables of the h-th feeder.

In summary, to apply the proposed method, matrices ${\mathit{U}}_{j,k}^h$ and ${\mathit{T}}_{j,k}^{h,i}$ in Equation (25) can be evaluated by the following procedure:

• Calculate the Jacobian matrices ${\mathbf{J}}^{sub}$ for the MV/LV supplying system and ${\mathbf{J}}_j^h$ for each j-th branch of each h-th feeder, from the analytic derivatives of Equation (1) evaluated in the initial operating point;
• Calculate the matrices ${\mathbf{M}}_j^h$ and ${\mathbf{N}}_{j,k}^h$ for each j-th node of each h-th feeder, according to Equations (5) and (6);
• Calculate the matrices ${\mathbf{A}}_k^h$ and the vectors ${\mathit{a}}^{\mathit{h}}$ through ${\overline{\mathbf{A}}}_k^h$, ${\overline{\mathit{a}}}^h$, which are provided by Equation (11), for each h-th feeder;
• Calculate the matrices ${\mathbf{B}}_{j,k}^h$ and the vectors ${\mathit{b}}_j^h$ for each j-th node of each h-th feeder according to Equation (14);
• Calculate the vector $\mathbf{\alpha }$ through $\overline{\mathbf{\alpha }}$ which is provided by Equation (23);
• Calculate the matrices ${\mathit{U}}_{j,k}^h$ and ${\mathit{T}}_{j,k}^{h,i}$ according to Equation (26), for each j-th node of each h-th feeder.

Finally, by substituting Equation (22) into (15) the variable $\Delta {\mathit{x}}_0^h$ can be expressed in the same form as Equation (25); then, the MV/LV substation variables can be obtained by using Equations (19) and (20).

3. Case Study

Reference is made to the 0.4 kV distribution system in Figure 5, supplied by an MV/LV substation with a 20/0.4 kV transformer of 0.25 MVA rated power. The 20 kV distribution system is represented by its Thevenin equivalent as seen from the MV/LV substation, assuming a 1000 MVA short-circuit power and an open-circuit voltage $V_{MV}$ = 1.0 per unit (p.u.). Concerning the MV/LV transformer, it is assumed that $X_{tr}$ = 0.06 p.u. and $R_{tr}$ = 0.005 p.u.. Two 0.4 kV radial feeders are connected to the LV substation busbar (m = 2). Both of the feeders are composed of $n_1$ = $n_2$ = 7 branches, which are characterized by the same impedances, reported in

Table 1. Electrical parameters and loads for the branches of both the feeders. p.u.: per unit.

From Node	To Node	Both Feeders		Feeder 1		Feeder 2
		R (p.u.)	X (p.u.)	PL (p.u.)	QL (p.u.)	PL (p.u.)	QL (p.u.)
0	1	0.0105	0.0025	0.0832	0.0416	0.0928	0.0464
1	2	0.0060	0.0014	0.8400	0.3780	0.2756	0.1376
2	3	0.0114	0.0027	0.0600	0.0300	0.0872	0.0436
3	4	0.0080	0.0011	0.0	0.0	0.0872	0.0436
4	5	0.0095	0.0014	0.1776	0.0800	0.2756	0.1376
5	6	0.0052	0.0007	0.0600	0.0300	0.2220	0.1112
6	7	0.0040	0.0006	0.1200	0.0600	0.0572	0.0284

. In the same Table 1, the active and reactive powers absorbed by the uncontrolled loads are also reported for each feeder. In the remainder, all p.u. quantities are referred to a 25 kVA power basis. Four DERs are connected to the LV distribution network, which can inject both active and reactive powers: ${\mathrm{DER}}_1$ and ${\mathrm{DER}}_2$ are connected to nodes 4 and 6 of feeder 1, respectively; ${\mathrm{DER}}_3$ and ${\mathrm{DER}}_4$ are connected to nodes 2 and 5 of feeder 2.

The case study is organized in four parts. Firstly, an application of the proposed method is explained step by step. Secondly, a model validation is performed by comparing the results given by the proposed model with the exact load flow solutions. Thirdly, the results of the proposed method are compared with a perturb-and-observe method in terms of both computational time and accuracy. Finally, an application of the proposed model is presented to quantify the impact of the active and reactive powers injected by DERs on nodal voltages.

3.1. Model Application

A load flow solution is evaluated by MATPOWER [19] (version 5.1, Power System Engineering Research Center (PSERC), Arizona State University, Tempe, AZ, USA) considering the distribution system in the initial operating condition in which all the power injections by the DERs are null. Starting from this initial point, the proposed method is applied and the sensitivity matrices ${\mathit{U}}_{j,k}^h$ and ${\mathit{T}}_{j,k}^{h,i}$ in Equation (25) are evaluated. Such sensitivity matrices allow for evaluating the impact of any injection by DERs on out-flowing powers and on voltage amplitudes of each node in the distribution system, as explained in the following for two examples.

Firstly, let the connection of only two DERs be considered, namely ${\mathrm{DER}}_1$ at node 4 of feeder 1 and ${\mathrm{DER}}_4$ at node 5 of feeder 2 (Figure 5). Referring to the shunt active and reactive powers at the DERs connection nodes, $\Delta P_{S,4}^1$ and $\Delta Q_{S,4}^1,$ respectively, indicate the active and reactive powers injected by ${\mathrm{DER}}_1$, whereas $\Delta P_{S,5}^2$ and $\Delta Q_{S,5}^2$, respectively, indicate the active and reactive powers injected by ${\mathrm{DER}}_4$. The impact of the two DERs on the electrical variables of the j-th node belonging to feeder 1 is derived from Equation (25):

$$\Delta {\mathit{x}}_j^1={\mathit{U}}_{j,4}^1\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right)+{\mathit{T}}_{j,5}^{1,2}\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right).$$

(27)

Matrix ${\mathit{U}}_{j,4}^1$ represents the linear dependence of the variations of the electrical variables of the nodes of feeder 1 on the power injection by ${\mathrm{DER}}_1$, which is installed along the same feeder 1, whereas matrix ${\mathit{T}}_{j,5}^{1,2}$ represents the linear dependence of the variations of the electrical variables of the nodes of feeder 1 on the power injection by ${\mathrm{DER}}_4$, which is installed along the other feeder 2. This latter dependence is due to the coupling between the feeders introduced by the MV/LV supplying system. Considering, for example, the third node of feeder 1 ($j=3$), the variations of the electrical variables are evaluated by:

$$\left(\begin{array}{c}\Delta P_3^1\\ \Delta Q_3^1\\ \Delta {V^2}_3^1\end{array}\right)=\left(\begin{array}{ccc}-1.0065& -0.0030& 0\\ -9.0·{10}^{-4}& -1.0004& 0\\ 0.0589& 0.0262& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right)+\left(\begin{array}{ccc}-5.0·{10}^{-6}& -4.5·{10}^{-5}& 0\\ -7.8·{10}^{-7}& -6.4·{10}^{-6}& 0\\ 0.0015& 0.0123& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right).$$

As expected, the values of the coefficients of ${\mathit{T}}_{3,5}^{1,2}$ show that ${\mathrm{DER}}_4$ has a negligible impact on the out-flowing powers but a significant impact on the voltage amplitude of node 3.

Similar considerations about the sensitivity matrices can be made when considering the variations of the electrical variables of any generic j-th node belonging to feeder 2 caused by the injections of the two DERs:

$$\Delta {\mathit{x}}_j^2={\mathit{U}}_{j,5}^2\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right)+{\mathit{T}}_{j,4}^{2,1}\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right).$$

(28)

Considering, for example, the third node of feeder 2 ($j=3$), the variations of the electrical variables are evaluated by:

$$\left(\begin{array}{c}\Delta P_3^2\\ \Delta Q_3^2\\ \Delta {V^2}_3^2\end{array}\right)=\left(\begin{array}{ccc}-1.0241& -0.0119& 0\\ -0.0034& -1.0017& 0\\ 0.0603& 0.0269& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right)+\left(\begin{array}{ccc}-1.5·{10}^{-5}& -1.3·{10}^{-4}& 0\\ -2.1·{10}^{-6}& -1.8·{10}^{-5}& 0\\ 0.0014& 0.0122& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right).$$

Secondly, let the connection of all four DERs be considered (Figure 5). Similarly to Equation (27), the impact of the four DERs on the electrical variables of the j-th node belonging to feeder 1 is derived from Equation (25):

$$\Delta {\mathit{x}}_j^1={\mathit{U}}_{j,4}^1\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right)+{\mathit{U}}_{j,6}^1\left(\begin{array}{c}\Delta P_{S,6}^1\\ \Delta Q_{S,6}^1\\ 0\end{array}\right)+{\mathit{T}}_{j,2}^{1,2}\left(\begin{array}{c}\Delta P_{S,2}^2\\ \Delta Q_{S,2}^2\\ 0\end{array}\right)+{\mathit{T}}_{j,5}^{1,2}\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right).$$

(29)

Considering, for example, the third node of feeder 1 ($j=3$), the variations of the electrical variables are evaluated by:

$$\left(\begin{array}{c}\Delta P_3^1\\ \Delta Q_3^1\\ \Delta {V^2}_3^1\end{array}\right)=\left(\begin{array}{ccc}-1.0065& -0.0030& 0\\ -9.0·{10}^{-4}& -1.0004& 0\\ 0.0589& 0.0262& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,4}^1\\ \Delta Q_{S,4}^1\\ 0\end{array}\right)+\left(\begin{array}{ccc}-1.0162& -0.0077& 0\\ -0.0023& -1.0011& 0\\ 0.05962& 0.0265& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,6}^1\\ \Delta Q_{S,6}^1\\ 0\end{array}\right)+\left(\begin{array}{ccc}-5.0·{10}^{-6}& -4.5·{10}^{-6}& 0\\ -7.0·{10}^{-7}& -6.4·{10}^{-6}& 0\\ 0.0013& 0.0122& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,2}^2\\ \Delta Q_{S,2}^2\\ 0\end{array}\right)+\left(\begin{array}{ccc}-5.5·{10}^{-6}& -4.5·{10}^{-5}& 0\\ -7.8·{10}^{-7}& -6.4·{10}^{-6}& 0\\ 0.0015& 0.0123& 0\end{array}\right)\left(\begin{array}{c}\Delta P_{S,5}^2\\ \Delta Q_{S,5}^2\\ 0\end{array}\right).$$

3.2. Model Validation

Starting from the initial operating condition in which all the power injections by the DERs are null, the proposed method is applied and the sensitivity matrices ${\mathit{U}}_{j,k}^h$ and ${\mathit{T}}_{j,k}^{h,i}$ in Equation (25) are evaluated. Then, new operating conditions are considered, characterized by active $P_{DER}$ and reactive $Q_{DER}$ power injections by the DERs. In particular, two different conditions are considered:

• Case A: $P_{DER,g}=10$ kW and $Q_{DER,g}=5$ kVAr for $g=1,\dots ,4$,
• Case B: $P_{DER,g}=20$ kW and $Q_{DER,g}=10$ kVAr for $g=1,\dots ,4$.

In both of the cases, the variations of nodal voltages and active and reactive power flows along the two feeders are evaluated by Equation (25) and added to the corresponding initial values, yielding the final values. These latter ones are compared with the results obtained by solving the load flow using MATPOWER in both Cases A and B. The results are compared in Figure 6 for feeder 1 (a) and 2 (b). In each figure, the first graph reports the nodal voltage amplitudes, and the second and the third graphs, respectively, the branch active and reactive power flows along the considered feeder. In each graph, the results obtained by the proposed method and by the load-flow solution are plotted for Cases A and B.

Considering feeder 1 (Figure 6a), in Case A, the evaluated voltage profile presents a decreasing trend from the substation LV busbar until node 2, and, then, it inverts its trend until node 6, in which, again, the voltage presents a slight decrease towards the end of the feeder. The reason for such a voltage profile is to be found in the power flows reported in the other graphs, which present two inversions, one between nodes 1 and 2, and the other one between nodes 5 and 6. In Case B, the voltage profile is strictly increasing along feeder 1 except for the last branch, because of the negative values of active and reactive powers, that is power flows toward the supplying substation. It is worth noting that the powers flowing out at the last nodes of the feeders are always null because the proposed model imposes the border constraints (8). In both of the cases and in all of the graphs, it is apparent that the results obtained by the proposed method accurately follow the same trend as the ones obtained by the load-flow solution. Similar considerations can be made analyzing feeder 2 (Figure 6b). Since loading conditions and DER allocations differ with respect to feeder 1, different voltage profiles and power flows are obtained. In Case A, power flows always remain positive, and, consequently, the voltage profile is strictly decreasing along the feeder. In Case B, power flows are negative from the LV substation until node 4 and then positive from node 5 until the end of the feeder; consequently, the voltage profile is increasing until node 5 and then slightly decreasing. In addition, for this feeder, the results obtained by the proposed method present the same trends as the ones obtained by load-flow solution.

To quantify the accuracy of the results with respect to the ones obtained by load-flow solution, the maximum relative errors on the variations of the nodal voltage amplitudes and on the branch active and reactive power flows along each feeder are reported in

Table 2. Maximum relative errors.

Case	Feeder Number	Maximum Errors (%) on the Variations of
		Voltage Amplitudes	Active Power Flows	Reactive Power Flows
A	1	3.3	4.4	1.9
	2	3.1	2.8	1.3
B	1	6.5	8.5	3.6
	2	5.5	5.4	2.4

for Cases A and B. The relative errors are always limited to between 1 and 9 %. The largest errors are obtained in Case B and feeder 1, which actually represents an extreme operating condition as clearly evidenced by the voltage profile in Figure 6.

3.3. Comparison with a Perturb-and-Observe Method

The proposed method is compared with an alternative method that is used for sensitivity analysis, namely a perturb-and-observe method of the type in [15]. In this approach, the values of the system variables in an assigned initial operating condition of the distribution system are evaluated by solving a load flow. Then, a small DER injection, firstly of active power and secondly of reactive power, for one node at a time, is added to the initial operating condition and a new load-flow is solved. Finally, the sensitivity coefficients of the system variables with respect to each DER injection are evaluated as the ratio between the variations of the system variables with respect to their corresponding values in the initial operating condition and the amplitude of the considered DER injection. In particular, for the considered distribution system, the same initial operating condition as the one adopted for the proposed method is assumed, in which all the power injections by the DERs are null. Then, small injections of active power (1 kW) and of reactive power ($1$kVAr) are imposed on one node at a time, for all the 14 nodes along the feeders: the variations of the system variables and the related sensitivity coefficients are obtained by 28 load-flow solutions.

Comparing the proposed method with the perturb-and-observe one in terms of computation times of the sensitivity coefficients, both methods require the solution of a load-flow in the same initial operating condition. The computation time of this load-flow solution is assumed as reference time unit and indicated as LF. The proposed method evaluates the coefficient matrices in Equation (25) using the closed-form expressions derived in Section 2 and the related computation time is equal to 0.20 LF. On the other hand, the perturb-and-observe method evaluates the sensitivity coefficients by 28 load-flow solutions with a computational time equal to 17.6 LF (it is less than 28 LF because of computational efficiency improvements that can be introduced). In conclusion, even for a small network as the considered one, the proposed method requires a computational time smaller of two orders of magnitude than the one required by a perturb-and-observe method.

To compare the proposed method with the perturb-and-observe one in terms of accuracy of the results, Case B defined in Section 3.2 is considered. Firstly, the voltage variations due to the DER injections are evaluated by using, respectively, the sensitivity coefficients for the perturb-and-observe method and Equation (25) for the proposed method. Then, the voltage variations are added to the corresponding values of the nodal voltages in the initial operating condition, yielding the new voltage amplitudes that account for the DER injections. In Figure 7, such voltage amplitudes evaluated by the proposed method, and by the perturb-and-observe method, are reported for both feeders and compared with the exact values provided by a load-flow solution in Case B. From the comparison, it is apparent that the results provided by the proposed method are more accurate than the ones obtained by the perturb-and-observe method. The main reason is that the proposed method uses the analytical derivative functions rather than the numerical derivatives used by the perturb-and-observe method.

3.4. Distributed Energy Resources Impact Evaluation

In the following, an example of the application of the proposed method to the evaluation of the impact of DERs on the nodal voltages of the distribution system is described to give evidence of the effectiveness of the method. The injection by ${\mathrm{DER}}_2$ connected to node 6 of feeder 1 and by ${\mathrm{DER}}_3$ connected to node 2 of feeder 2 are considered. The advantage of the proposed method is its linearity that allows for separately analyzing the impact in terms of nodal voltage variations of active and reactive power injection by each DER, to separately analyze the impact of the injections $P_{DE{R}_2}$, $Q_{DE{R}_2}$, $P_{DE{R}_3}$, $Q_{DE{R}_3}$. Then, the impact of any combination of these injections can be obtained by simply adding the voltage variations that each considered injection produces at each node, provided that all voltage variation is referred to the same initial operating condition, which is always assumed to be the one in which all the power injections by DERs are null.

Referring to active power injection, Figure 8a reports the variations of the nodal voltage amplitudes along both feeders caused by imposing $P_{DE{R}_2}=$ 10 kW and evaluated by the proposed method. The injection has a different impact on the two feeders. In particular, on feeder 1, to which ${\mathrm{DER}}_2$ is connected, the power injection causes a significant variation of the nodal voltage profile, which increases starting from the LV busbar until node 6, and then remains quite constant until the end of the feeder. Concerning feeder 2, the power injection along the other feeder has a practically null effect on the nodal voltage profile. In conclusion, the active power injection has a large impact on the feeder to which the DER is connected; on the contrary, it has a negligible impact on other feeders because of its negligible impact on the voltage drop along the MV/LV transformer.

Referring to reactive power injection, Figure 8b reports the variations of the nodal voltage amplitudes along both feeders caused by imposing $Q_{DE{R}_2}=$ 10 kVAr and evaluated by the proposed method. Concerning feeder 1, the reactive power injection causes a significant variation of the nodal voltage profile, which increases from the LV busbar until node 6, and then remains quite constant until the end of the feeder. Concerning feeder 2, there is a significant variation on all the nodal voltages, but the shape of the voltage profile does not change, since all the nodes are subject to approximately the same voltage variation. In conclusion, the reactive power injection has a significant impact on the feeder to which the DER is connected, but also on other feeders, because it causes a reduction of the voltage drop along the MV/LV transformer.

Comparing the graphs (a) and (b) in Figure 8, it is apparent that $P_{DE{R}_2}$ has a higher impact on the nodal voltages of feeder 1 with respect to $Q_{DE{R}_2}$ because the branch resistances $R_j$ are larger than the branch reactances $X_j$ (Table 1); on the contrary, $Q_{DE{R}_2}$ has a higher impact on the nodal voltages of feeder 2 with respect to $P_{DE{R}_2}$ because the substation reactance $X_{eq}$ is larger than the substation resistance $R_{tr}$.

Similar considerations can be made for Figure 9, which reports the voltage variations caused by imposing $P_{DE{R}_3}=10$ kW (a) and $Q_{DE{R}_3}=10$ kVAr (b). In this case, the major impact is on the voltages of feeder 2, where ${\mathrm{DER}}_3$ is connected.

Comparing Figure 8 and Figure 9, it is apparent that the impact on the voltage profile of ${\mathrm{DER}}_2$ on feeder 1 is higher than the impact of ${\mathrm{DER}}_3$ on feeder 2 because the former ${\mathrm{DER}}_2$ is connected near the end of feeder 1, whereas the latter ${\mathrm{DER}}_3$ is connected near the beginning of feeder 2. On the contrary, the impacts of ${\mathrm{DER}}_2$ on feeder 2 and of ${\mathrm{DER}}_3$ on feeder 1 are very similar, which is negligible for active power injections and about $0.0024$ p.u. for reactive power injections. Then, the impact of a DER on the voltage profile of a feeder, other than the one at which the DER itself is connected, does not depend on the position of the connection node.

Using the voltage variations that quantify the impact of each DER injection, it is possible to:

• Identify the VMs with respect to voltage regulation problems: e.g., from the graphs (a) in Figure 8 and Figure 9, it is apparent that the two feeders belong to two different VMs as far as active power injections are considered;
• Quantify the effectiveness of each DER on the voltage regulation in a VM: e.g., from the comparison of Figure 8 and Figure 9, it is apparent that ${\mathrm{DER}}_2$ has a major effect on the voltage regulation of a VM along feeder 1 with respect to ${\mathrm{DER}}_3$;
• Quantify the interaction between the voltage regulation devices that are installed in a VM [17]: e.g., if a VM is defined along feeder 1, the voltage variation at node 4 caused by ${\mathrm{DER}}_2$ and the voltage variation at node 6 caused by ${\mathrm{DER}}_1$ represent the interaction between the voltage regulators acting on ${\mathrm{DER}}_1$ and ${\mathrm{DER}}_2$.

4. Conclusions

A new model of LV distribution systems with DERs has been proposed, which is suitable for aggregating an existing distribution system into a cluster of virtual microgrids. The model is based on a linearized formulation of the branch flow equations and of the MV/LV supplying system. Then, the model is solved by imposing the border constraints at the MV busbar of the substation and at the end-nodes of the feeders and the coupling equations at LV substation busbar. The resulting closed form solution expresses the variations of the electrical variables of the LV distribution network as linear functions of the DERs injections. A numerical case study has evidenced the satisfactory accuracy of the results obtained by the proposed method with respect to the exact load flow solution. A comparison of the proposed method with a perturb-and-observe method is discussed with reference to computation time and accuracy of the results, thus giving evidence of the better performance of the proposed method. Eventually, an application to the evaluation of the impact of DERs on the voltage regulation problem of an LV distribution system has been presented, quantifying the sensitivity of the nodal voltages to the active and reactive power injections. Future work will extend the model to account for more complex network topologies and for unbalanced operating condition.