Distributed reactive power feedback control for voltage regulation and loss minimization

I Introduction

Recent technological advances, together with environmental and economic challenges, have been motivating the deployment of small power generators in the low voltage and medium voltage power distribution grid. The availability of a large number of these generators in the distribution grid can yield relevant benefits to the network operation, which go beyond the availability of clean, inexpensive electrical power. They can be used to provide a number of ancillary services that are of great interest for the management of the grid .

We focus in particular on the problem of optimal reactive power compensation for power losses minimization and voltage regulation. In order to properly command the operation of these devices, the distribution network operator is required to solve an optimal reactive power flow (ORPF) problem. Powerful solvers have been designed for the ORPF problem, and advanced optimization techniques have been recently specialized for this task . However, this approach assumes that an accurate model of the grid is available, that all the grid buses are monitored, that loads announce their demand profiles in advance, and that generators and actuators can be dispatched on a day-ahead, hour-ahead, and real-time basis. For this reason, these solvers are in general offline and centralized, and they collect all the necessary field data, compute the optimal configuration, and dispatch the reactive power production at the generators.

These tools cannot be applied directly to the ORPF problem faced in low/medium voltage power distribution networks. The main reasons are that not all the buses of the grid are monitored, individual loads are unlikely to announce they demand profile in advance, the availability of small size generators is hard to predict (being often correlated with the availability of renewable energy sources). Moreover, the grid parameters, and sometimes even the topology of the grid, are only partially known, and generators are expected to connect and disconnect, requiring an automatic reconfiguration of the grid control infrastructure (the so called plug and play approach).

Different strategies have been recently proposed in order to address these issues. Purely local algorithms have been proposed, in which each generator is operated according to its own measurements , in order to compensate for the voltage rise caused by its own active power injection . Because of the absence of coordination between microgenerators, the full potential of the microgenerators for voltage regulation is not exploited in these strategies .

Different coordination strategies have been then proposed, for example by casting the problem into the framework of resource allocation and by using hierarchical dispatch schemes . A two-stage approach has been proposed in , where microgenerators first attempt to regulate their voltage autonomously, and they involve their neighbors in this task if their regulation capability is not sufficient.

Finally, some distributed approaches that do not require any central controller, but still require measurements at all the buses of the distribution grid, have been proposed. In order to derive a distributed algorithm for this problem, different convex relaxation methods have been applied, and various distributed optimization algorithms have been specialized for the resulting convex ORPF problem .

Only recently, algorithms that are truly scalable in the number of generators and do not require the monitoring of all the buses of the grid, have been proposed for the problem of power loss minimization (with no voltage constraints) . While these algorithms have been designed by specializing classical nonlinear optimization algorithms to the ORPF problem, they can also be considered as feedback control strategies. Indeed, the key feature of these algorithms is that they require the alternation of measurement and actuation based on the measured data, and therefore they are inherently online algorithms. In particular, the reactive power injection of the generators is adjusted by these algorithms based on the phasorial voltage measurements that are performed at the buses where the generators are connected. The resulting closed loop system features a tight dynamic interconnection of the physical layer (the grid, the generators, the loads) with the cyber layer (where communication, computation, and decision happen). In this paper, we design a distributed feedback algorithm for the ORPF problem with voltage constraints, in which the goal is to minimize reactive power flows while ensuring that the voltage magnitude across the network lies inside a given interval, and that the microgenerators’ reactive power limits are not violated. The analysis of the convergence of this strategy is based on an assumption of homogeneity of the X/R ratio of the power lines across the network. The robustness of the proposed solution with respect to possible variability in these parameters has been investigated via simulations.

In Section III, a cyber-physical model for a smart power distribution grid is provided. In Section V, the ORPF problem with voltage constraints is formulated. A feedback control strategy for its solution is derived in Section VI, by using the tools of dual decomposition. A synchronous and an asynchronous version of the algorithm are presented in Section VI and Section VII, respectively. The convergence of both the proposed algorithms is studied in Section VIII. Some simulations are provided in Section IX, while Section X concludes the paper discussing some relevant features of the feedback nature of the proposed strategy.

II Mathematical preliminaries and notation

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\sigma,\tau)$ be a directed graph, where $\mathcal{V}$ is the set of nodes, $\mathcal{E}$ is the set of edges, and $\sigma,\tau:\mathcal{E}\rightarrow\mathcal{V}$ are two functions such that edge $e\in\mathcal{E}$ goes from the source node $\sigma(e)$ to the terminal node $\tau(e)$.

Given two nodes $h,k\in\mathcal{V}$, we define the path $\mathcal{P}_{hk}$ as the sequence of adjacent nodes, without repetitions, that connect node $h$ to node $k$.

Let $A\in\{0,\pm 1\}^{{|\mathcal{E}|}\times n}$ be the incidence matrix of the graph $\mathcal{G}$, defined via its elements

If the graph $\mathcal{G}$ is connected (i.e. for every pair of nodes there is a path connecting them), then ${\mathbf{1}}$ is the only vector in the null space $\ker A$, ${\mathbf{1}}$ being the column vector of all ones. We define by ${\mathbf{1}}_{v}$ the vector whose value is $1$ in position $v$, and everywhere else.

III Cyber-physical model of a smart power distribution grid

In this work, we envision a smart power distribution network as a cyber-physical system, in which

the physical layer consists of the power distribution infrastructure, including power lines, loads, microgenerators, and the point of connection to the transmission grid;

the cyber layer consists of intelligent agents, dispersed in the grid, and provided with actuation, sensing, communication, and computational capabilities.

For the purpose of this paper, we model the physical layer of a smart power distribution network as a directed graph $\mathcal{G}$, in which edges represent the power lines, and nodes represent buses (see Figure 1, middle panel). Buses correspond to loads, microgenerators, and also the point of connection to the transmission grid (called point of common coupling, or PCC, and indexed as node ).

We limit our study to the steady state behavior of the system, when all voltages and currents are sinusoidal signals at the same frequency. Each signal can therefore be represented via a complex number $y=|y|e^{j\angle y}$ whose absolute value $|y|$ corresponds to the signal root-mean-square value, and whose phase $\angle y$ corresponds to the phase of the signal with respect to an arbitrary global reference. In this notation, the steady state of the grid is described by the following system variables (see Figure 1, lower panel):

We model the grid power lines as series impedances, neglecting their shunt admittance. For every edge $e$ of the graph, we define by $z_{e}$ the impedance of the corresponding power line. We assume the following.

All the power lines in the grid have the same inductance/resistance (X/R) ratio, but possibly different impedance magnitude, i.e.

for any $e$ in $\mathcal{E}$ and for a fixed $\theta$.

This assumption is satisfied when the X/R ratio of the power lines of the grid is relatively homogeneous, which is reasonable in many practical cases (see for example the IEEE standard testbeds ). In Section IX we investigate what is the effect of a possible variability of the X/R ratio, and we show how the proposed strategy is very robust against this possible uncertainty.

Under this assumption, we have a linear relation between bus voltages and currents in the form

where $L:=A^{T}Z^{-1}A$ is the weighted Laplacian of the graph, in which $A$ is the incidence matrix of $\mathcal{G}$, and $Z=\operatorname{diag}(|z_{e}|,e\in\mathcal{E})$ is the diagonal matrix of the magnitudes of line impedances.

Each node $v$ of the grid is then characterized by a law relating its injected current $i_{v}$ with its voltage $u_{v}$. We model the PCC as an ideal sinusoidal voltage generator at the microgrid nominal voltage $U_{N}$ with arbitrary fixed angle $\psi$

In the power system analysis terminology, node is then a slack bus with fixed voltage magnitude and angle.

We model loads and microgenerators (that is, every node $v$ of the microgrid except the PCC) via the following law relating the voltage $u_{v}$ and the current $i_{v}$

where $s_{v}$ is the injected complex power. The quantities

are denoted as active and reactive power, respectively. The complex powers $s_{v}$ corresponding to grid loads are such that $\{p_{v}<0\}$, meaning that positive active power is supplied to the devices. The complex powers corresponding to microgenerators, on the other hand, are such that $\{p_{v}\geq 0\}$, as positive active power is injected into the grid. In the power system analysis terminology, all nodes but the PCC are being modeled as constant power or P-Q buses. Microgenerators fit in this model, as they generally are commanded via a complex power reference and they can inject it independently from the voltage at their point of connection .

III-B Cyber layer

We assume that every microgenerator, and also the PCC, correspond to an agent in the cyber layer (upper panel of Figure 1). We denote by ${\mathcal{C}}$ (with $|{\mathcal{C}}|=m$) this subset of the nodes of $\mathcal{G}$. Each agent is provided with some computational capability, and with some sensing capability, in the form of a phasor measurement unit (i.e. a sensor that can measure voltage amplitude and angle ). Agents that corresponds to a microgenerator can also actuate the system, by commanding a set point for the amount of reactive power injected by that microgenerator (see Figure 2).

Finally, agents can communicate, via some communication channel that could possibly be the same power lines (via power line communication). Motivated by this possibility, we define the neighbors in the cyber layer in the following way.

Let $h\in{\mathcal{C}}$ be an agent of the cyber layer. The set of agents that are neighbors of $h$, denoted as $\mathcal{N}(h)$, is the subset of ${\mathcal{C}}$ defined as

We assume that every agent $h\in{\mathcal{C}}$ knows its set of neighbors $\mathcal{N}(h)$, and can communicate with them. Notice that this architecture can be constructed by each agent in a distributed way, for example by exploiting the PLC channel (as suggested for example in ). This allows also a plug-and-play reconfiguration of such architecture when new agents are connected to the grid.

IV Approximated model and its properties

In this section we review an approximate explicit solution of the nonlinear equations (1), (2), and (3) which has been proposed in . This approximation will play a crucial role in deriving our distributed control strategy to solve the optimal reactive power flow problem that we will introduce in next section. In order to present the approximated solution, we need the following technical lemma.

The matrix $X$ depends only on the topology of the grid power lines and on their impedance. The matrix $X$ has some notable properties, including the fact that

where $Z^{\text{eff}}_{hk}$ represents the effective impedance of the power lines between node $h$ and $k$. Notice that, if the grid is radial (i.e. $\mathcal{G}$ is a tree) then $Z^{\text{eff}}_{hk}$ is simply the impedance of the only path from node $h$ to node $k$. Now let us introduce the following block decomposition of the vector of voltages $u$

By adopting this block decomposition as before, we have

Consider the physical model described by the set of nonlinear equations (1), (2), and (3). Node voltages then satisfy

where the little-o notation means that $\lim_{U_{N}\rightarrow\infty}\frac{o(f(U_{N}))}{f(U_{N})}=0$.

It descends directly from Proposition 1 in . ∎

The quality of this approximation relies on having large nominal voltage $U_{N}$ and relatively small currents injected by the inverters (or supplied to the loads). This assumption is verified in practice, and corresponds to correct design and operation of power distribution networks, where indeed the nominal voltage is chosen sufficiently large (subject to other functional constraints) in order to deliver electric power to the loads with relatively small power losses on the power lines. The model proposed in Proposition 1 extends the DC power flow model [28, Chapter 3] to the case in which lines are not purely inductive. This way, the model is able to describe the voltage drop on the lines, and therefore also the corresponding power losses, in a form that is conveniently linear in the complex power injections and demands. We conclude this section by introducing the following matrix $G$.

The following matrix $G$ satisfies the conditions.

The proof of uniqueness, that we omit here, follows exactly the same steps as in the proof of Lemma 4. ∎

The matrix $G$ has also a remarkable sparsity pattern, as the following lemma states.

The matrix $G$ has the sparsity pattern induced by the Definition 1 of neighbor agents in the cyber layer, i.e.

The proof is provided in the Appendix A, where we also discuss how the elements of $G$ can be estimated by the agents, given a local knowledge of the power grid parameters.

V Optimal reactive power flow problem

We consider the problem of commanding the reactive power injection of the microgenerators in order to minimize power distribution losses on the power lines and to guarantee that the voltage magnitude and the reactive power injection stay within pre-assigned intervals. The decision variables (or, equivalently, the inputs of the system) are therefore the reactive power setpoints $q_{h},h\in{\mathcal{C}}\backslash\{0\}$, compactly written as $q_{G}$.

Power distribution losses can be expressed as a function of the voltage drop on the lines and therefore, in a matricial quadratic form, as

Given a lower bound $U_{\text{min}}$ and an upper bound $U_{\text{max}}$ for the voltage magnitudes, and a lower bound $q_{\text{min}}$ and an upper bound $q_{\text{max}}$ for the reactive power injected by each microgenerator, we can therefore formulate the following optimization problem,

where voltages $u$ are a function of the decision variables $q_{G},$ via the implicit relation defined by the system of nonlinear equations (1), (2), and (3). From a control design prospective, the system-wide problem that we are considering is therefore characterized by

the measured output variables $\begin{bmatrix}u_{0}\\ u_{G}\end{bmatrix}$,

the unmeasured disturbances $p_{L}$, $q_{L}$, $p_{G}$.

The goal of this paper is to design a control algorithm to tackle the ORPF problem in a distributed fashion, where each microgenerator $h$ is allowed to communicate only with its neighbors in the cyber layer, i.e., the agents in $\mathcal{N}(h)$.

While the decision variables of the ORPF problem (i.e. the input variables $q_{G}$) do not include the reactive power supplied by the PCC (i.e. $q_{0}=u_{0}\bar{i}_{0}$), this quantity will also change every time the reactive power setpoints of the generators are updated by the algorithm, because the inherent physical behavior of the slack bus (the PCC) ensures that equations (1), (2) and (3) are satisfied at every time.

In the above formulation we have assumed that all the microgenerators have to satisfy the same reactive power injection constraint. This scenario can be seamlessly extended to the case of heterogeneous microgenerators, where $q_{\text{h,min}}\leq q_{h}\leq q_{\text{h,max}}$ being the values $q_{\text{h,min}}$, $q_{\text{h,max}}$, $h\in{\mathcal{C}}$, in general different for different microgenerators.

VI A synchronous algorithm based on dual decomposition

In this section, in order to design a distributed feedback control strategy to solve the ORPF problem, we apply the tool of dual decomposition to (11). Specifically, we use the approximate explicit solution of the nonlinear equations (1), (2), and (3) introduced in Proposition 1, to derive update steps for a dual ascent algorithm that can be implemented distributively by the agents and that can be used as a feedback control update law.

It is convenient to reformulate the problem (11) via a change of coordinates, obtaining

Basically, in (12f) with the respect to (11f), we have squared and normalized the constraints on the voltage magnitude and we have normalized the constraints on the power injection. While these modifications does not have any effect on the optimization problem, they will allow us to simplify the derivation of the algorithm we are going to present. The Lagrangian of the problem (12) is

where $\lambda_{\text{min}}$, $\lambda_{\text{max}}$, $\mu_{\text{min}}$, $\mu_{\text{max}}$ are the Lagrangian multipliers (i.e. the dual variables of the problem) and $u$, $v_{G}$, $w_{G}$ are functions of the decision variables $q_{G}$, even if the dependence has been omitted. To have a more compact notation let

A dual ascent algorithm consists in the iterative execution of the following alternated steps

dual gradient ascent step on the dual variables

where the $[\cdot]_{+}$ operator corresponds to the projection on the positive orthant, and where $\gamma$ is a suitable positive constant;

minimization of the Lagrangian with respect to the primal variables $q_{G}$

Observe that the updates of the Lagrange multipliers can be performed naturally in a distributed way by the agents, based on their local measurement of the violation of the voltage and power constraints. Indeed let $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, $\mu_{\text{min},h}$ and $\mu_{\text{max},h}$ be the components of the the Lagrange multipliers $\lambda_{\text{min}}$, $\lambda_{\text{max}}$, $\mu_{\text{min}}$ and $\mu_{\text{max}}$, respectively, related to the compensator $h$. Then it easily follows that the dual step can be implemented as

The crucial point is to derive an expression for the minimizer $q_{G}(t+1)$ in (14) that can be computed distributively by the compensators. To do so, we exploit the approximation introduced in Proposition 1 obtaining a value of $q_{G}(t+1)$ which is equivalent to the one in (14) up to a term which vanishes to zero for large nominal voltage $U_{N}$.

The resulting update corresponds to the algorithm that we now present. We assume here that the agents are coordinated, i.e., they can update their state variables $q_{h}$, $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, $\mu_{\text{min},h}$ and $\mu_{\text{max},h}$, synchronously.

Let all agents (except the PCC) store the auxiliary scalar variables $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, $\mu_{\text{min},h}$ and $\mu_{\text{max},h}$. Let $\gamma$ be a positive scalar parameter, and let $\theta$ be the impedance angle defined in Assumption 1. Let $G_{hk}$ be the elements of the sparse matrix $G$ defined in Lemma 2. At every synchronous iteration of the algorithm, each agent $h\in{\mathcal{C}}\backslash\{0\}$ executes the following operations in order:

it measures its voltage $u_{h}$ and it gathers the voltage measurements

it updates the auxiliary variables $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, $\mu_{\text{min},h}$, $\mu_{\text{max},h}$, as

it gathers from its neighbors the updated values of the Lagrange multipliers $\mu_{\text{min},k}$, $\mu_{\text{max},k}$, $k\in\mathcal{N}(h)$;

based on the new values of $\lambda_{\text{max},h}$, $\lambda_{\text{min},h}$ and of $\mu_{\text{min},k}$, $\mu_{\text{max},k}$, $k\in\mathcal{N}(h)$, it updates the injected reactive power $q_{h}$ as

Observe that the above algorithm can be implemented in a completely distributed fashion. Indeed each agent is required to exchange information only with its neighbors in the cyber layer.

The following Proposition shows how the update (15) approximates the primal step (14).

Consider the synchronous algorithm above described. Then

namely, the update (15), minimizes the Lagrangian with respect to the primal variables, up to a term that vanishes for large $U_{N}$.

The details of the proof of Proposition 2 are postponed to Appendix B. It is based on the following technical lemma, which will be useful again later.

Consider the Lagrangian $\mathcal{L}(q_{G},\lambda)$ defined in (13). The partial derivative with respect to the primal variables $q_{G}$ is

Also the proof of Lemma 4 is in Appendix B.

Via these steps, we therefore specialized the dual ascent steps to the ORPF problem that we are considering, and we obtained a distributed feedback control law for the system. We study the convergence of the closed loop system in Section VIII.

It is important to notice that the proposed synchronous algorithm requires that the agents actuate the system at every iteration, by updating the set point for the amount of reactive power injected by the microgenerators. Only by doing so, the subsequent measurement of the voltages will be informative of the new state of the system. The resulting control strategy is thus a feedback strategy that necessarily requires the real-time interaction of the controller (the cyber layer) with the plant (the physical layer), as depicted in Figure 4. This tight interaction between the cyber layer and the physical layer is the fundamental feature of the proposed approach, and allows to drive the system towards the optimal configuration, which in principle depends on the reactive power demands of the loads, without collecting this information from them. In a sense, the algorithm is inferring this hidden information from the measurement performed on the system during its execution.

VII Asynchronous algorithm

In order to avoid the burden of system-wide coordination among the agents, we also propose an asynchronous version of the algorithm, in which the agents corresponding to the microgenerators update their state $(q_{h},\lambda_{\text{max},h},\lambda_{\text{min},h},\mu_{\text{max},h},\mu_{\text{min},h})$ independently one from the other, based on the information that they can gather from their neighbors.

We assume that each agent (except for the agent located at the PCC) is provided with an individual timer, by which it is triggered, and that no coordination is present between these timers: they tick randomly, with exponentially, identically distributed waiting times.

Let all agents (except the PCC) store four auxiliary scalar variables $\lambda_{\text{max},h},\lambda_{\text{min},h},\mu_{\text{max},h},\mu_{\text{min},h}$. Let $\gamma$ be a positive scalar parameter, and let $\theta$ be the impedance angle defined in Assumption 1. Let $G_{hk}$ be the elements of the matrix $G$ defined in Lemma 2.

When agent $h\in{\mathcal{C}}\backslash\{0\}$ is triggered by its own timer, it performs the following actions in order:

it measures its voltage $u_{h}$ and it gathers from its neighbors the voltage measurements

and the values of the Lagrange multipliers

it updates the auxiliary variables $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, $\mu_{\text{min},h}$, $\mu_{\text{max},h}$, as

based on the new value of $\lambda_{\text{min},h}$, $\lambda_{\text{max},h}$, it updates the injected reactive power $q_{h}$ as

The update equations for the asynchronous algorithm are exactly the same of the synchronous case. Here, however, we update both the primal and the dual variable of the agents independently and asynchronously. Also the analysis of the convergence of this algorithm is postponed to the next section.

VIII Convergence analysis

In this section, we investigate the convergence of both the synchronous algorithm proposed in Section VI and of the asynchronous algorithm proposed in Section VII.

In order to do so, we rewrite the terms that appeared in the dual ascent update step, namely

using the expression introduced in Proposition 1 for the voltages. We start from

By plugging in the approximate solution (8), via some algebraic manipulations, we can express $v_{G}$ as

The proposed dual ascent step can therefore be rewritten in compact form as

The update step for the primal variables can be rewritten based on Proposition 2 and Lemma 4, obtaining

In the analysis that follows, we study the approximated description of the closed loop system in which we neglect the infinitesimal terms. Notice that, by doing so, both the voltages $u$ and the squared voltage magnitudes $v_{G}$ become affine functions of the decision variables $q_{G}$. By plugging those expressions in the formulation of the ORPF problem (12), one obtains the following strictly convex quadratic problem with linear inequality constraints

for which strong duality holds. The rest of the section is split into two subsections: in the first one we consider the synchronous version of the algorithm, in the second one the asynchronous version.

For the synchronous version of the algorithm, we consider the update equations

for the primal variables. Observe that (23) and (24) differ from (19) and from (21) only by infinitesimal terms, and they correspond to the standard equation for the dual ascent steps for (22). Indeed, the equilibrium $(q_{G}^{*},\nu^{*})$ of (23)-(24) is characterized by

which correspond to the necessary conditions for the optimality according to Uzawa’s saddle point theorem .

It will be useful in the following to define $\sigma_{\text{min}}$ and $\sigma_{\text{max}}$ as the minimum and the maximum eigenvalue of $M$, respectively. The following result characterizes the convergence of the algorithm described by (23) and (24).

Consider the optimization problem (22) and the dynamic system described by the update equations (23) and (24). Then the trajectory $t\to q(t)$ converges to the optimal primal solution $q_{G}^{*}$ if

We conclude this subsection by specializing the above result to the case where either only voltage constraints or only power constraints are considered. Observe that if we take into account only voltage constraints then the matrix $\Phi$ and the vector $b$ become

and only the multipliers $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ are employed in the algorithm, while if we consider only power constraints then

and only the multipliers $\mu_{\text{min}}$ and $\mu_{\text{max}}$ are needed. The following results follow from Theorem 1.

Consider the optimization problem (22), where $\Phi$ and $b$ are given as in (25), and the dynamic system described by the update equations (23) and (24). Then the trajectory $t\to q(t)$ converges to the optimal primal solution $q_{G}^{*}$ if

Consider the optimization problem (22), where $\Phi$ and $b$ are given as in (26), and the dynamic system described by the update equations (23) and (24). Then the trajectory $t\to q(t)$ converges to the optimal primal solution $q_{G}^{*}$ if

VIII-B Asynchronous case

Let us define the random sequence $h(t)\in{\mathcal{C}}\backslash\{0\}$ which tells which agent has been triggered at iteration $t$ of the algorithm. Because of Assumption 2, the random process $h(t)$ is an i.i.d. uniform process on the alphabet ${\mathcal{C}}\backslash\{0\}$. If we repeat the same analysis, neglecting the infinitesimal terms, we obtain the following update equations for the primal and dual variables, instead of (23) and (24). In these equations, only the component $h(t)$ of the vectors $\lambda_{\text{min}}$, $\lambda_{\text{max}}$, $\mu_{\text{min}}$, $\mu_{\text{max}}$ and $q_{G}$ is updated at time $t$, namely,

Notice that, also in the asynchronous case, Uzawa’s necessary conditions for optimality are satisfied at the equilibrium of (27), (28) and (29). For the asynchronous version of the algorithm we can provide theoretical results only when $\Phi$ and $b$ assume the form in (25) or in (26), namely, when we consider either only voltage constraints or only power constraints. However in the numerical section we show the effectiveness of the asynchronous algorithm when $\Phi$ and $b$ assume the general form in (20).

Consider the optimization problem (22), where $\Phi$ and $b$ are given as in (25), and the dynamic system described by the update equations (27) and (28) (for the multipliers $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$) and (29). Let Assumption 2 hold. Then the evolution $t\to q(t)$ converges almost surely to the optimal primal solution $q_{G}^{*}$ if

Consider the optimization problem (22), where $\Phi$ and $b$ are given as in (26), and the dynamic system described by the update equations (27) and (28) (for the multipliers $\mu_{\text{min}}$ and $\mu_{\text{max}}$) and (29). Let Assumption 2 hold. Then the evolution $t\to q(t)$ converges almost surely to the optimal primal solution $q_{G}^{*}$ if

IX Simulations

The algorithm has been tested on the testbed IEEE 37 , which is an actual portion of 4.8kV power distribution network located in California. The load buses are a blend of constant-power, constant-current, and constant-impedance loads, with a total power demand of almost 2 MW of active power and 1 MVAR of reactive power (see for the testbed data). The length of the power lines range from a minimum of 25 meters to a maximum of almost 600 meters. The impedance of the power lines differs from edge to edge (for example, resistance ranges from 0.182 $\Omega$/km to 1.305 $\Omega$/km). However, the inductance/resistance ratio exhibits a smaller variation, ranging from $X/R=0.5$ to $0.67$. This justifies Assumption 1, in which we claimed that $\angle z_{e}$ can be considered constant across the network. We considered the scenario in which $5$ microgenerators have been deployed in this portion of the power distribution grid (see Figure 5).

The lower bound for voltage magnitudes has been set to 4700 V. Both the synchronous and the asynchronous algorithm presented in Section VI and VII have been simulated on a nonlinear exact solver of the grid . The approximate model presented in Proposition 1 has not been used in these simulations, being only a tool for the design of the algorithm and for the study of the algorithm’s convergence.

A time-varying profile for the loads has been generated, in order to simulate the effect of slowly varying loads (e.g. the aggregate demand of a residential neighborhood), fast changing / intermittent demands (e.g. some industrial loads).

The results of the simulation have been plotted in Figure 6 for the asynchronous case, while the synchronous case has not been reported, being very similar. In order to tune the parameter $\gamma$, based on the similarity between the conditions in Propositions 3 and 4, and those in Corollaries 1 and 2, we conjecture that the bound derived in Theorem 1 is also valid in the asynchronous case. We have therefore chosen $\gamma$ to be one half of such bound. The power distribution losses, the lowest voltage magnitude measured by the microgenerators, and the reactive power injection of one of the microgenerators, are reported. The dashed line represents the case in which no reactive power compensation is performed. The thick black line represents the best possible strategy that solves the ORPF problem (11) (computed via a numerical centralized solver that have real time access to all the grid parameters and load data). The thin red line represents the behavior of the proposed algorithm.

It can be seen that the proposed algorithm achieves practically the same performance of the centralized solver, in terms of power distribution losses. Notice however that the proposed algorithm does not have access to the demands of the loads, which are unmonitored. The agents, located only at the microgenerators, can only access their voltage measurements and share them with their neighbors. Notice moreover that, as expected for duality based methods, the voltage constraints can be momentarily violated. Therefore, in the time varying case simulated in this example, the voltage sometimes falls slightly below the prescribed threshold, when the power demand of the loads present abrupt changes. It should be remarked, however, that the extent of this constraint violation depends on the rate at which the algorithm is executed, compared with the rate of variation of loads, and on the fact that an exact (and thus aggressive) primal update step has been implemented. The same behavior cannot be observed for the power constraints, as the reactive power set-point has been saturated in order to simulate the typical implementation of power inverters, which cannot accept set-point references that exceed their rated power. Notice that the reactive power reference is almost constant when voltage constraints are not active (as the primal step is exact, and therefore the algorithm reaches the optimal point immediately). When the constraints are active, the evolution depends instead on the update of the Lagrange multipliers.

Finally, in Figure 7, we investigated the robustness of the algorithm with respect to possible larger variations of the X/R ratio of the power lines. The same testbed has been modified in order to have X/R ratios ranging from 0.36 to 2.6. Despite the fact that Assumption 1 is needed for the technical results of the paper, simulations show how the effect on the closed loop behavior of the controlled systems is minimal. Intuitively, this is due to the feedback nature of the control strategy: as the violation of the constraints is integrated in the feedback loop (see Figure 4), this violation is guaranteed to go to zero, as long as the closed loop system is stable.

X Conclusions

In this paper we proposed a distributed control law for optimal reactive power flow in a smart power distribution grid, based on a feedback strategy. Such a strategy requires the interleaving of actuation and sensing, and therefore the control action (the reactive power injections $q_{h},h\in{\mathcal{C}}\backslash\{0\}$) is a function of the real time measurements (the voltages $u_{h},h\in{\mathcal{C}}$). According to this interpretation, the active power injections in the grid ($p_{h},h\in\mathcal{V}$) and the reactive power injection of the loads ($q_{h},h\in\mathcal{V}\backslash{\mathcal{C}}$) can be considered as disturbances for the control system. As explained in Remark 3, these quantities do not need to be known to the controller, and the agents are implicitly inferring them from the measurements. It is also well known that the presence of feedback in the control action makes the closed loop behavior of the system less sensitive to model uncertainties, as shown in the simulations. These features differentiate the proposed algorithm from most of the ORPF algorithms available in the power system literature, with the exception of some works, like , where however the feedback is only local, with no communication between the agents, and of and . Moreover, in the proposed feedback strategy, the controller does not need to solve any model of the grid in order to find the optimal solution. The computational effort required for the execution of the proposed algorithm is therefore minimal. These features are extremely interesting for the scenario of power distribution networks, where real time measurement of the loads is usually not available, and the grid parameters are partially unknown.

While a feedback approach to the ORPF problem is a recent approach, similar methodologies have been used to solve other tasks in the operation of power grids (see Figure 8). In particular, in order to achieve realtime power balance of demand and supply, synchronous generators are generally provided with a local feedback control that adjusts the input mechanical power according to frequency deviation measurements (the primary frequency control in , see also ). By adding a communication channel (a cyber layer) that enables coordination among the agents, it is possible to drive the system to the configuration of minimum generation costs . Notice that, in this scenario, generators do not have access to the aggregate active power demand of the loads, but infer it from the purely local frequency measurements. In this sense, this example share some qualitative similarities with the original approach presented in this paper.

As suggested in , a control-theoretic approach to optimization problems (including ORPF) enables a number of analyses on the performance of the closed loop system that are generally overlooked. Examples are $L_{2}$-like metrics for the resulting losses in a time-varying scenario (e.g. the preliminary results in ), robustness to measurement noise and parametric uncertainty, stability margin against communication delays. These analyses, still not investigated, are also of interest for the design of the cyber architecture, because they can provide specifications for the communication channels, communication protocols, and computational resources that need to be deployed in a smart distribution grid.

Appendix A g-parameters

Let us define the following parameters $g_{hk}$, $h,k\in{\mathcal{C}}$.

For each pair $h,k\in{\mathcal{C}}$, let us define the parameter

i.e. the current that would be injected at node $h$ if

node $k$ was replaced with a unitary voltage generator;

all other nodes in ${\mathcal{C}}$ (the other agents) were replaced by short circuits;

all nodes not in ${\mathcal{C}}$ (load buses) were replaced by open circuits.

Notice that the parameters $g_{hk}$ depend only on the grid electric topology, and that

Figure 9 gives a representation of this definition. Notice that, in the special case in which the paths from $h$ to its neighbors are all disjoint and unique paths, then $G_{hk}=(\sum_{e\in\mathcal{P}_{hk}}|z_{e}|)^{-1}$, i.e. the inverse of the impedance of the path connecting $h$ to $k$.

As suggested in , these parameters can be estimated in an initialization phase via some ranging technologies over the PLC channel. Alternatively, this limited amount of knowledge of the grid topology can be stored in the agents at the deployment time. Finally, the same kind of information can be also inferred by specializing the procedures that use the extended capabilities of the generator power inverters for online grid sensing and impedance estimation .

The following Lemma shows that the elements $g_{hk}$ in Definition 2 correspond to the elements $G_{hk}$ of the matrix $G$ in Lemma 2.

Let $G$ be the matrix defined in Lemma 2. Then, for all $h,k$, $G_{hk}=g_{hk}$ as defined in Definition 2.

Let $\hat{G}$ be the matrix whose elements are the parameters $g_{hk}$. From Definition 2, we have that when $i_{L}=0$,

From circuit theory considerations, this implies that $\hat{G}{\mathbf{1}}=0$. From (1) and by using the matrix $X$ defined in Lemma 4, if $i_{L}=0$, we have

By plugging (31) into (32), we obtain $\left(I-{\mathbf{1}}{\mathbf{1}}_{0}^{T}\right)=\left[\begin{smallmatrix}0&0\\ 0&M\end{smallmatrix}\right]\hat{G}$. ∎

Appendix B Proof of Lemma 4 and Proposition 2

By using the fact that $\bar{u}^{T}Lu=u^{\prime T}Lu^{\prime}+u^{\prime\prime T}Lu^{\prime\prime}$, we have

where we used the fact that $L{\mathbf{1}}=0$ and that, by Lemma 4, $LX=I-{\mathbf{1}}_{0}{\mathbf{1}}^{T}$.

The same approximate solution (8), via some algebraic manipulations, allows us to express $v_{G}$ as

and finally, from (33), (34), (36) and (37),

It can be shown, by using Lemma 2 and via some algebraic manipulation, that the update (15) can be also rewritten as

which, by using the expression for $u$ provided by Proposition 1, is equal to

Then, after the update, by plugging the former into the expression for the partial derivative of the Lagrangian with the respect to $q_{G}$, provided in Lemma 4, we obtain

and therefore the update minimized the Lagrangian with respect to the primal variables, up to a term that vanishes for large $U_{N}$. ∎

Appendix C Proof of Theorem 1, Corollaries 1 and 2, and Propositions 3 and 4

It is straightforward to see that the dual of (22) is

Since (22) is quadratic optimization problem that we have assumed feasible and the constraint is expressed by a linear affine inequality, the Slater’s condition [41, p. 226] holds and then there is zero duality gap between (22) and (39). Observe that, by plugging (24) into (23) we obtain

that is the update of $\nu$ is a projected gradient ascent algorithm for the dual function $g(\lambda)$. Then any optimal solution $\nu^{*}$ of (39) is a fixed point for (41) and satisfies

while the primal optimal solution, from (24), has the form

It is worth to notice that (40) has not necessarily a unique solution: given a particular solution $\nu_{1}^{*}$, if there exists $v\in\ker(\Phi^{T})$ such that $\nu^{*}_{2}=\nu_{1}^{*}+v=[\nu^{*}_{2}]_{+}$, then also $\nu^{*}_{2}$ is an optimal solution. Despite that, $q^{*}_{G}$ is unique. In fact we have

Notice that we have, $\forall\nu_{1},\nu_{2}\geq 0$, that

and then the gradient $\nabla g(\nu)$ is Lipschitz continuous with Lipschitz constant equals to $\|\Phi M^{-1}\Phi^{T}\|$. Then, from Prop. 2.3.2 in , if

the algorithm (41) converges to a maximizer of $g(\nu)$ and then we reach the optimal solution $q_{G}^{*}$ of the primal optimization problem. We have

Being $\Phi M^{-1}\Phi^{T}$ a positive semi-definite symmetric matrix, its norm is equal to its spectral radius. It can be shown that the spectrum of $\Phi M^{-1}\Phi^{T}$ is given by $\Lambda(\Phi M^{-1}\Phi^{T})=\{0\}\cup\Lambda(2\Xi)$ where

The characteristic polynomial of $\Xi$ is

where $\sigma_{j}$ is the $j$-th eigenvalues of $M$ and where, without loss of generality we assume that $0<\sigma_{1}\leq\sigma_{2}\leq\ldots\leq\sigma_{m-1}$. Thus we can see that

Consider the case where only voltage constraints are considered. Then we have that

from which it follows that $\rho\left(\Phi M^{-1}\Phi^{T}\right)=2\sin^{2}\theta\sigma_{\text{max}}$.

Instead when only power constraints are taken into account we have that

from which we get that $\rho\left(\Phi M^{-1}\Phi^{T}\right)=2\sigma_{\text{min}}^{-1}$. ∎

Consider the update equations (27), (28) and (29) for the dual variables $\nu$ and for the primal variables $q_{G}$. Let $(q_{G}^{*},\nu^{*})$ be a solution of the optimization problem, which satisfies (42) and the KKT conditions

We introduce the following two quantities

Without loss of generality let us assume that node $h$ is the node performing the update at the $t$-th iteration. The update for the variable $x$ is given by

where we used (29) and (46a). Now, let us consider first the case where only voltage constraints are taken into account. Via some algebraic manipulations we can write from (27) that

Thanks to the fact that $|a_{+}-b_{+}|\leq|a-b|$ we can write that

Now observe that Assumption 2 implies there exists almost surely a positive integer $T$ such that any node has performed an update within the window $[0,T]$. Moreover observe from (C) that $x_{h}(t+1)=-{\mathbf{1}}_{h}^{T}[-\sin\theta I\,\,\sin\theta I]^{T}y(t+1)$. It follows that, for $t\geq T$, $x(t)=-M^{-1}\Phi^{T}y(t)$. Hence we can write

Let us adopt the decomposition $y=y_{\perp}+y_{\parallel}$ where $y_{\perp}\perp\ker\Phi^{T}$ and $y_{\parallel}\in\ker\Phi^{T}$. It follows

If $c=0$, then $y$ goes to which implies that also $x$ goes to zero and then $q_{G}$ tends to the optimal primal solution $q_{G}^{*}$. Otherwise if $c\neq 0$, we have from (49) that

This implies that the trajectory $\nu(t)$ tends to the set

This implies the boundedness of the sequence $\nu(t)$ and the convergence of $q_{G}(t)$ to $q_{G}^{*}$, being from (29)

We now briefly repeat similar steps for the case where only power constraints are taken into account. In this case we have from (27) that

From the expression for $x_{h}$ in (C), it follows that

Recalling that, in this case, $\Phi=[-I\,\,\,I]^{T}$, from the above inequalities we get that, as in (48), that

Again the convergence of $q_{G}(t)$ to the optimal solution $q_{G}^{*}$ is guaranteed if condition (50) is satisfied. ∎

References