Active Distribution Grids offering Ancillary Services in Islanded and Grid-connected Mode

TL;DR

This paper explores how active distribution grids with distributed energy resources can provide frequency control services and operate reliably in islanded mode, ensuring grid stability and resilience.

Contribution

It introduces a centralized optimal power flow model that enables ADGs to offer ancillary services while maintaining islanded operation capabilities.

Findings

The scheme can respond to various frequency control products.

It guarantees energy capacity for critical loads in islanded mode.

The approach is feasible on low-voltage distribution networks.

Abstract

Future active distribution grids (ADGs) will incorporate a plethora of Distributed Generators (DGs) and other Distributed Energy Resources (DERs), allowing them to provide ancillary services in grid-connected mode and, if necessary, operate in an islanded mode to increase reliability and resilience. In this paper, we investigate the ability of an ADG to provide frequency control (FC) in grid-connected mode and ensure reliable islanded operation for a pre-specified time period. First, we formulate the operation of the grid participating in European-type FC markets as a centralized multi-period optimal power flow problem with a rolling horizon of 24 hours. Then, we include constraints to the grid-connected operational problem to guarantee the ability to switch to islanded operation at every time instant. Finally, we explore the technical and economic feasibility of offering these services on a balanced low-voltage distribution network. The results show that the proposed scheme is able to offer and respond to different FC products, while ensuring that there is adequate energy capacity at every time step to satisfy critical load in the islanded mode.

Equations

$$\underset{\mathbf{u}}{\min }\sum _{t=t_{\text{MPC}}}^{t_{\text{MPC}}+T}\left(C_{\text{t}}^{\text{curt,g}}+C_{\text{t}}^{\text{curt,l}}+C_{\text{t}}^{\text{exc}}+C_{\text{t}}^{\text{AS}}\right)\Delta t,$$

$$C_{\text{t}}^{\text{curt,g}}=\sum _{j\in \mathcal{R}}{c}_{\text{t}}^{\text{curt,g}}\cdot (P_{\text{j,t}}^{\text{curt,g}}+f_{\text{isl}}\cdot P_{\text{j,t}}^{\text{curt,g,isl}}),$$

$$C_{\text{t}}^{\text{curt,l}}=\sum _{j\in \mathcal{L}}{c}_{\text{t}}^{\text{curt,l,isl}}\cdot P_{\text{j,t}}^{\text{curt,l,isl}},$$

$$C_{\text{t}}^{\text{exc}}=c_{\text{t}}^{\text{buy}}\cdot P_{\text{1,t}}^{\text{buy}}-c_{\text{t}}^{\text{sell}}\cdot P_{\text{1,t}}^{\text{sell}},$$

$$C_{\text{t}}^{\text{AS}}=c_{\text{t}}^{\text{bid}}\cdot P_{\text{bid}},$$

$$P_{\text{j,t}}^{\text{inj}}=P_{\text{j,t}}^{\text{g}}-P_{\text{j,t}}^{\text{lflex}}-\left(P_{\text{j,t}}^{\text{B,ch}}-P_{\text{j,t}}^{\text{B,dis}}\right),$$

$$Q_{\text{j,t}}^{\text{inj}}=Q_{\text{j,t}}^{\text{g}}+Q_{\text{j,t}}^{\text{B}}-P_{\text{j,t}}^{\text{lflex}}\cdot \tan ({\varphi }_{\text{l}}),$$

$$I_{\text{j,t}}^{\text{inj}}=\left(\frac{(P_{\text{j,t}}^{\text{inj}}+j{Q}_{\text{j,t}}^{\text{inj}}{)}^{\ast }}{{\bar{V}}_{\text{j,t}}^{\ast }}\right),$$

$$I_{\text{t}}^{\text{br}}=BIBC\cdot I_{\text{t}}^{\text{inj}},$$

$$\Delta V_{\text{t}}=BCBV\cdot I_{\text{t}}^{\text{br}},$$

$$V_{\text{j,t}}=V_{\text{slack}}-\Delta V_{\text{tap}}\cdot {\rho }_{\text{t}}+\Delta V_{\text{t}},$$

$${\rho }_{min}\le {\rho }_{\text{t}}\le {\rho }_{max},$$

$$|I_{\text{i,t}}^{\text{br}}|\le I_{\text{i}}^{\text{max}},$$

$$V_{\text{min}}\le |V_{\text{j,t}}|\le V_{\text{max}}.$$

$$V_{\text{min}}$$

$$P_{\text{j,t}}^{\text{g,min}}\le P_{\text{j,t}}^{\text{g}}\le P_{\text{j,t}}^{\text{g,max}},Q_{\text{j,t}}^{\text{g,min}}\le Q_{\text{j,t}}^{\text{g}}\le Q_{\text{j,t}}^{\text{g,max}},$$

$$P_{\text{j,t}}^{\text{lflex}}=P_{\text{j,t}}^{\text{l}}+f_{\text{j,t}}^{\text{lflex}}\cdot P_{\text{j,t}}^{\text{shift}},$$

$$-1\le f_{\text{j,t}}^{\text{lflex}}\le 1,$$

$$\sum _{t=t_{\text{0}}}^{t_{\text{MPC}}-1}{f}_{\text{j,t}}^{\text{lflex}}+\sum _{t=t_{\text{MPC}}}^{t_{\text{MPC}}+T}{f}_{\text{j,t}}^{\text{lflex}}=0,$$

$$So{C}_{\text{min}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}\le E_{\text{j,t}}^{\text{B}}\le So{C}_{\text{max}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}},$$

$$E_{\text{j,1}}^{\text{B}}=E_{\text{j,t0}},$$

$$E_{\text{j,t}}^{\text{B}}=E_{\text{j,t-1}}^{\text{B}}+({\eta }_{\text{B}}\cdot P_{\text{j,t}}^{\text{B,ch}}-\frac{P_{\text{j,t}}^{\text{B,dis}}}{{\eta }_{\text{B}}})\cdot \Delta t,$$

$$0\le P_{\text{j,t}}^{\text{B,ch}}\le P_{\text{j,max}}^{\text{B}},0\le P_{\text{j,t}}^{\text{B,dis}}\le P_{\text{j,max}}^{\text{B}},$$

$$P_{\text{j,t}}^{\text{B,ch}}+P_{\text{j,t}}^{\text{B,dis}}\le \text{max}(P_{\text{j,t}}^{\text{B,ch}},P_{\text{j,t}}^{\text{B,dis}}),$$

$$|Q_{\text{j,t}}^{\text{B}}|\le \max \left\{P_{\text{j,t}}^{\text{B,ch}},P_{\text{j,t}}^{\text{B,dis}}\right\}\cdot tan({\varphi }_{\text{max}}^{\text{B}}),$$

$$\sum _{j\in \mathcal{B}}\left(P_{\text{max,j}}^{\text{B}}-P_{\text{j,t}}^{\text{B,dis}}+P_{\text{j,t}}^{\text{B,ch}}\right)\ge P_{\text{bid}},$$

$$\sum _{j\in \mathcal{B}}\left(P_{\text{max,j}}^{\text{B}}-P_{\text{j,t}}^{\text{B,ch}}+P_{\text{j,t}}^{\text{B,dis}}\right)\ge P_{\text{bid}},$$

$$\sum _{j\in \mathcal{B}}\left(E_{\text{j,t}}^{\text{B}}-So{C}_{\text{min}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}\right)\ge P_{\text{bid}}\cdot \Delta t_{\text{1}},$$

$$\sum _{j\in \mathcal{B}}\left(So{C}_{\text{max}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}-E_{\text{j,t}}^{\text{B}}\right)\ge P_{\text{bid}}\cdot \Delta t_{\text{1}},$$

$$\sum _{j\in \mathcal{B}}\left(P_{\text{max,j}}^{\text{B}}-P_{\text{j,t}}^{\text{B,dis}}+P_{\text{j,t}}^{\text{B,ch}}\right)+\sum _{j\in \mathcal{R}}\left(P_{\text{j,t}}^{\text{g,max}}-P_{\text{j,t}}^{\text{g}}\right)\ge P_{\text{bid}},$$

$$\sum _{j\in \mathcal{B}}\left(P_{\text{max,j}}^{\text{B}}-P_{\text{j,t}}^{\text{B,ch}}+P_{\text{j,t}}^{\text{B,dis}}\right)+\sum _{j\in \mathcal{R}}\left(P_{\text{j,t}}^{\text{g}}\right)\ge P_{\text{bid}},$$

$$\begin{gathered} E_{\text{j,t+ϑ}}^{\text{B,2,+}}=\underset{\text{shortstackPrevious BESS energy content}}{\underbrace{E_{\text{j,t+ϑ-1}}^{\text{B,2,+}}}}-\underset{\text{shortstackIf all energy had to be provided by the BESS}}{\underbrace{\frac{1}{{\eta }_{\text{B}}}\cdot P_{\text{bid}}\cdot \Delta t_{\text{2,ϑ}}^{\text{+}}}}+ \\ +\underset{\text{shortstackPart that can be provided by PVs (max. 50\% of worst case call)}}{\underbrace{min\left\{0.5\cdot \frac{1}{{\eta }_{\text{B}}}\cdot P_{\text{bid}},\sum _{\mathcal{R}}\left(P_{\text{j,t+ϑ}}^{\text{g,max}}-P_{\text{j,t+ϑ}}^{\text{g}}\right)\right\}\cdot \Delta t_{\text{2,ϑ}}^{\text{+}}}}+ \\ +\underset{\text{shortstackScheduled BESS operation}}{\underbrace{{\eta }_{\text{B}}\cdot P_{\text{j,t+ϑ}}^{\text{B,ch}}\cdot \Delta t-\frac{1}{{\eta }_{\text{B}}}\cdot P_{\text{j,t+ϑ}}^{\text{B,dis}}\cdot \Delta t}} \end{gathered}$$

$$\begin{gathered} E_{\text{j,t+ϑ}}^{\text{B,2,-}}=E_{\text{j,t+ϑ-1}}^{\text{B,2,-}}+{\eta }_{\text{B}}\cdot P_{\text{bid}}\cdot \Delta t_{\text{2,ϑ}}^{\text{-}} \\ -min\left\{0.5\cdot {\eta }_{\text{B}}\cdot P_{\text{bid}},\sum _{\mathcal{R}}{P}_{\text{j,t+ϑ}}^{\text{g}}\right\}\cdot \Delta t_{\text{2,ϑ}}^{\text{-}} \\ +{\eta }_{\text{B}}\cdot P_{\text{j,t+ϑ}}^{\text{B,ch}}\cdot \Delta t-\frac{1}{{\eta }_{\text{B}}}\cdot P_{\text{j,t+ϑ}}^{\text{B,dis}}\cdot \Delta t \end{gathered}$$

$$So{C}_{\text{min}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}\le E_{\text{j,t+ϑ}}^{\text{B,2,+}}\le So{C}_{\text{max}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}},$$

$$So{C}_{\text{min}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}\le E_{\text{j,t+ϑ}}^{\text{B,2,-}}\le So{C}_{\text{max}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}.$$

$$\begin{gathered} \sum _{\text{j}\in \mathcal{B}}\left(P_{\text{j,max}}^{\text{B, j}}-P_{\text{j,t}}^{\text{B,dis}}+P_{\text{j,t}}^{\text{B,ch}}\right)+\sum _{\text{j}\in \mathcal{R}}\left(P_{\text{j,t}}^{\text{g,max}}-P_{\text{j,t}}^{\text{g}}\right)+ \\ \sum _{\text{j}\in \mathcal{F}}\left(f_{\text{j,t}}^{\text{lflex}}-(-1)\right)\cdot P_{\text{j,t}}^{\text{shift}}\ge P_{\text{bid}}, \end{gathered}$$

$$\begin{gathered} E_{\text{j,t+ϑ}}^{\text{B,3,+}}=E_{\text{j,t+ϑ-1}}^{\text{B,3,+}}-\frac{P_{\text{bid}}\cdot \Delta t_{\text{3}}}{{\eta }_{\text{B}}}+ \\ +min\{0.8\cdot \frac{1}{{\eta }_{\text{B}}}\cdot P_{\text{bid}},\sum _{\text{j}\in \mathcal{R}}\left(P_{\text{j,t}}^{\text{g,max}}-P_{\text{j,t}}^{\text{g}}\right)+ \\ +\sum _{\text{j}\in \mathcal{F}}\left(f_{\text{j,t}}^{\text{lflex}}+1\right)\cdot P_{\text{j,t}}^{\text{shift}}\}\cdot \Delta t_{\text{3}}+\left({\eta }_{\text{B}}\cdot P_{\text{j,t}}^{\text{B,ch}}-\frac{P_{\text{j,t}}^{\text{B,dis}}}{{\eta }_{\text{B}}}\right)\cdot \Delta t, \end{gathered}$$

$$So{C}_{\text{min}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}\le E_{\text{j,t+ϑ}}^{\text{B,3,+}}\le So{C}_{\text{max}}^{\text{B}}\cdot E_{\text{cap,j}}^{\text{B}}.$$

Full text

I Introduction

While moving towards a low-carbon, sustainable electricity system, future Distribution Networks (DNs) are expected to host a large share of Distributed Generators (DGs) to satisfy the demand currently supplied by fossil-fuel and nuclear power plants. DGs, coordinated with other Distributed Energy Resources (DERs), such as electric vehicles, Battery Energy Storage Systems (BESSs) and Flexible Loads (FLs), consequently amplify the role of DNs, making them an important part in ensuring grid reliability and resilience [1], and enabling them to provide ancillary services to transmission voltage levels [2]. Thus, it is crucial to operate modern DNs actively, i.e. controlling DERs to ensure secure, reliable and cost-effective operation.

In this paper, we consider a centralized method with existing communication infrastructure, which is a valid assumption in modern DNs that do not cover large geographical areas [3].

I-A Related work

Operating active DNs using optimization has been widely explored in literature, e.g. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Here, we only review work concerned with the DN capability to operate off-grid and the provision of ancillary services offered in grid-connected mode.

I-A1 Islanded operation

In [5], a Monte Carlo (MC) approach is applied in the design stage to determine the required BESS size to reliably operate in this mode but without incorporating the possibility of the BESS offering ancillary services in the connected mode. Reference [7] on the other hand optimizes the microgrid operational costs in grid-connected mode as a master problem, while ensuring islanded capability for multiple hours as a subproblem. However, neither the provision of energy based ancillary services nor the incorporation of voltage control, which requires the consideration of a network model, are addressed. A model-predictive-control (MPC) scheme, including the dynamics of the system, is used in [6] to predict future voltage instabilities and adjust the reactive power generation accordingly. Here, the focus is only on keeping voltages close to nominal values in islanded mode, not offering other grid-connected or islanded services. Further, [13] examines the behavior of a real BESS offering frequency control reserves and supporting islanded operation. The described setup uses a dispatchable diesel generator in addition to intermittent renewable energy sources but no network constraints are considered and the response of the BESS is based on heuristics rather than on centralized optimization.

I-A2 Grid-connected operation

In grid-connected mode, the main objective is usually to operate the DN in the most cost-effective way. A detailed review of the state-of-the-art research in microgrids is presented in [16], where the authors review around $400$ works, covering the areas of microgrid economics, operation, control, protection, and communications. Reference [4] investigates the economic evaluation of grid-connected microgrids that participate in real-time markets but without considering islanded operation. Ref. [15] focuses also on the optimal scheduling of an active DN providing frequency regulation, load leveling and ramping services. Offering ancillary services by various DER technologies is explored in [8], while [9, 10] analyze the economic feasibility and the potential amount of reserve provision by distributed generation. However, they do not consider BESSs or include network constraints. The technical feasibility of providing ancillary services with multiple microgrids as a pool bidder is investigated in [11], while [12] takes the perspective of a Transmission System Operator (TSO), minimizing its own expenditures and evaluating the competitiveness of DGs. Finally, several DSOs are already providing actual frequency control products to the TSO. For instance [13] and [14] discuss the operational experience of a BESS offering primary frequency control in the European interconnected network. However neither of [9, 10, 15, 13, 14] considers the network modeling within the optimization.

The consideration of grid-connected and islanded mode of an active distribution grid in combination with the provision of ancillary service taking into account also grid constraints is, to the best of our knowledge, not considered in any previous work.

I-B Contributions

In this paper, we propose a centralized optimization approach to operate an Active Distribution Grid (ADG). We explicitly incorporate uncertainty into the formulation and consider the opportunity of offering ancillary services. The proposed methodology is based on a multi-period, Chance-Constrained Optimal Power Flow (CC-OPF) formulation, where various frequency control products are offered by DERs. We ensure that at any point the DN can operate in islanded mode for a limited time by including additional constraints in the centralized problem. In this way, the DN operation considers both the uncertainties from RES as well as the potential need for islanded operation [17].

Consequently, the main contributions of this paper are:

• •

A multi-period CC-OPF formulation that:

• –

allows an ADG to offer ancillary services in grid-connected mode, while being able to switch to islanded mode at any time, and

• –

considers RES uncertainty through a rolling horizon strategy.

• •

A case-study analyzing various frequency control (FC) products and the performance of the proposed method in offering these products.

The remainder of the paper is organized as follows: Section II presents the mathematical formulation of the deterministic OPF considering ancillary service provision as well as islanded operation, while Section III accounts for uncertainty and presents the final CC-OPF Then, Section IV introduces the case study and simulation results for the islanded and grid-connected case. Finally, conclusions are drawn in Section V.

II Centralized deterministic OPF

In this section, we present the deterministic centralized OPF scheme used to compute the optimal DER setpoints. The objective considers both the grid-connected and the islanded mode simultaneously and is optimized in a rolling horizon fashion; we model FC products offered in the grid-connected mode as constraints, while at the same time enabling a potential switch to islanded mode for the following 24 hours.

II-A Centralized OPF

II-A1 Preliminaries

We consider a radial balanced distribution grid with $\mathcal{N}$ being the set of nodes using the index $j$, $\mathcal{T}$ the set of branches using the index $i$, $\mathcal{B}\subseteq\mathcal{N}$ the subset of nodes with BESS, $\mathcal{L}\subseteq\mathcal{N}$ the subset of nodes with loads, $\mathcal{F}\subseteq\mathcal{L}\subseteq\mathcal{N}$ the subset of flexible (controllable) loads, and $\mathcal{R}\subseteq\mathcal{N}$ the subset of nodes with DGs. The DER control measures (detailed below) are represented by the variable $u$, and the variables referring to the islanded operation mode by the subscript “isl”.

II-A2 Objective function

The objective function is defined as

[TABLE]

where $t_{\textrm{MPC}}$ denotes the current time step of the MPC algorithm, $T$ the rolling horizon period and $\Delta t$ the length of a time interval within the horizon.

At each discrete time $t$, the objective function consists of four terms;

$C_{\textrm{t}}^{\textrm{curt,g}}$

This term corresponds to the cost of generation curtailment in both the grid-connected and islanded mode and is given by

[TABLE]

where $P_{\textrm{j,t}}^{\textrm{curt,g}}=P_{\textrm{j,t}}^{\textrm{g,max}}-P_{\textrm{j,t}}^{\textrm{g}}$ is the curtailed power of the DG connected at node $j$ at time $t$ (resp., $P_{\textrm{j,t}}^{\textrm{curt,g,isl}}$ in the islanded case), $P_{\textrm{j,t}}^{\textrm{g,max}}$ the maximum available active power, and $P_{\textrm{j,t}}^{\textrm{g}}$ the actual in-feed; $c_{\textrm{t}}^{\textrm{curt,g}}$ is the cost of curtailment at time $t$, and $f_{\textrm{isl}}$ a constant scalar that adjusts the cost in the islanded case. The cost of generation curtailment is policy-related in the grid-connected case and can be very different from country to country. Typically, generators are compensated at the prevailing electricity market price, whereas in some European countries, they are compensated only for a small part of the curtailed energy [18]. In California, compensation for curtailment begins after a contractually agreed number of hours which vary among contracts [19]. However, in the islanded case, the operation of the DGs becomes more important, since they are the only sources to satisfy the local demand, i.e. no external grid is available. Thus, there is another value associated with the injection of power from DGs in the islanded case, which is accounted for by the scalar $f_{\textrm{isl}}$.

$C_{\textrm{t}}^{\textrm{curt,l}}$

This term represents the cost of load curtailment in the islanded mode and is given as

[TABLE]

where $c_{\textrm{t}}^{\textrm{curt,l,isl}}$ is the cost of load curtailment at time $t$, and $P_{\textrm{j,t}}^{\textrm{curt,l,isl}}$ the curtailed load. In the grid-connected case, we assume that any local generation-load mismatch can be covered from the transmission network without the need for load shedding in the grid-connected case, similar to [20].

$C_{\textrm{t}}^{\textrm{exc}}$

The third term includes the cost of exchanging power with the upper voltage levels and is given by

[TABLE]

where $c_{\textrm{t}}^{\textrm{buy}}$ ($c_{\textrm{t}}^{\textrm{sell}}$) is the price of buying (selling) electric energy from (to) the main grid. By considering different prices for buying and selling, i.e. at each time step buying electricity is more expensive than selling, we prioritize storing excess energy locally (promoting the self-consumption of the DN), over exporting power to higher voltage levels; $P_{\textrm{1,t}}^{\textrm{g}}=P_{\textrm{1,t}}^{\textrm{buy}}-P_{\textrm{1,t}}^{\textrm{sell}}$ ($P_{\textrm{1,t}}^{\textrm{buy}},P_{\textrm{1,t}}^{\textrm{sell}}\geq 0$) is the active power exchange measured at the substation making sure that the ADG cannot buy and sell electricity at the same time. A similar formulation is followed in [21] to determine the position, i.e. short or long, of an aggregator participating in energy markets.

$C_{\textrm{t}}^{\textrm{AS}}$

The final term corresponds to revenues from offering ancillary services to upper voltage levels, given by

[TABLE]

where $c_{\textrm{t}}^{\textrm{bid}}$ is the pay-as bid volume-weighted average price of the accepted bids in the frequency control market from the respective week of the previous year (assumed known) and $P_{\textrm{bid}}$ the bid (to be determined by the CC-OPF).

II-A3 Power balance constraints

The power injections at every node $j$ and time step $t$ are given by

[TABLE]

where $P_{\textrm{j,t}}^{\textrm{g}}$ and $Q_{\textrm{j,t}}^{\textrm{g}}$ are the active and reactive power injections of the generators at node $j$; $P_{\textrm{j,t}}^{\textrm{lflex}}$ and $P_{\textrm{j,t}}^{\textrm{lflex}}\cdot\tan(\phi_{\textrm{l}})$ are the active and reactive node demands (after control), with $cos(\phi_{\textrm{l}})$ being the power factor of the load; $Q_{\textrm{j,t}}^{\textrm{B}}$ is the reactive power of the BESS and, $P_{\textrm{j,t}}^{\textrm{B,ch}}$ and $P_{\textrm{j,t}}^{\textrm{B,dis}}$ are respectively the charging and discharging BESS active powers.

II-A4 Power flow constraints

In this work, we integrate the Backward/Forward Sweep (BFS) method into our power flow formulation [22, 23, 24]. The solution of the BFS power flow problem is achieved iteratively, by ”sweeping” the distribution network and updating the network variables at each iteration, which consists of two sweeps. First, in the backward sweep step, the current injections at all buses and the corresponding branch currents are calculated. Then, in the forward sweep step, the currents are used to calculate the voltage drop over all branches, updating the bus voltages for the next iteration of the algorithm. Within an OPF framework, we consider only one iteration to model network flows and to avoid the non-linearities introduced by the AC power flow equations. Subsequently, if the derived solution is not AC feasible, we update the voltages by projecting the solution into the AC feasible manifold [23], and re-run the OPF problem. This reformulation provides a sufficiently accurate approximation of the full AC OPF [25], is computationally tractable [24], and results in AC feasible solutions which can account for uncertainties (see [23] for more details). A single iteration of the BFS is used to replace the AC power-flow constraints in the OPF formulation as follows:

[TABLE]

where $\bar{V}_{\textrm{j,t}}^{*}$ is the voltage magnitude at node $j$ at time $t$,$~{}^{*}$ indicates the complex conjugate and the bar indicates that the value from the previous iteration is used; $I_{\textrm{t}}^{\textrm{inj}}=[I_{\textrm{j,t}}^{\textrm{inj}},j\in\mathcal{N}]$ and $I_{\textrm{t}}^{\textrm{br}}=[I_{\textrm{i,t}}^{\textrm{br}},i\in\mathcal{T}]$ represent respectively the vectors of bus injection and branch flow currents; $I_{\textrm{i,t}}^{\textrm{br}}$ is the $i$-th branch current; $BIBC$ is a matrix with ones and zeros, capturing the radial topology of the DN; the entries in $\Delta{V}_{\textrm{t}}$ correspond to the voltage drops over all branches and phases; $BCBV$ is a matrix with the complex impedances of the lines as elements; $V_{\textrm{slack}}$ is the voltage in per unit at the slack bus (here assumed to be $1\phase{0^{\circ}}$); $\Delta V_{tap}$ is the voltage magnitude change caused by one tap action of the On-Load Tap Changer (OLTC) transformer and assumed constant for all taps for simplicity; and, $\rho_{\textrm{t}}$ is an integer value defining the position of the OLTC position. The parameters ($\rho_{\textrm{min}},\rho_{\textrm{max}}$) are respectively the minimum and maximum tap positions of the OLTC transformer.

II-A5 Thermal loading and voltage constraints

The constraints for the current magnitudes at time $t$ are given by

[TABLE]

where $I_{\textrm{i}}^{\textrm{max}}$ is the maximum thermal limit for the $i$-th branch.

Similarly, the voltage constraints are given by

[TABLE]

where $(V_{\textrm{max}}$, $V_{\textrm{min}})$ are respectively the upper and lower acceptable voltage limits. However, the lower voltage magnitude limit results in a non-convex constraint [24]. By exploiting the fact that the voltage angles are typically small in distribution grids, we can approximate the complex voltage with its real part for the lower bound, as explained in [24]:

[TABLE]

II-A6 DER constraints

DG limits

In this work, we consider inverter-based DGs such as PVs. Their limits are given by

[TABLE]

where $P_{\textrm{j,t}}^{\textrm{g,min}}$, $P_{\textrm{j,t}}^{\textrm{g,max}}$, $Q_{\textrm{j,t}}^{\textrm{g,min}}$ and $Q_{\textrm{j,t}}^{\textrm{g,max}}$ are the upper and lower limits for active and reactive DG power at each node $j\in\mathcal{N}$ and time $t$. These limits vary depending on the type of the DG and the control schemes implemented. Usually, small DGs have technical or regulatory [26] limitations on the power factor they can operate at or reactive power they can produce. Any of these limitations can be captured in this constraint.

Controllable loads

We further consider flexible loads which can shift a limited amount of energy consumption in time. The loads are therefore modeled by

[TABLE]

where $P_{\textrm{j,t}}^{\textrm{shift}}$ is the shiftable load of the non-shiftable demand $P_{\textrm{j,t}}^{\textrm{l}}$; $f_{\textrm{j,t}}^{\textrm{lflex}}$ is the normalized factor defining the final load shift. The past values for $f_{\textrm{j,t}}^{\textrm{lflex}}$, i.e. for time instances $t=t_{\textrm{0}}$ (start of the simulation) to $t_{\textrm{MPC}}-1$, are constant. This is necessary due to the moving horizon approach and the fact that the total load at the end of the simulation period needs to be maintained which is ensured by (11c). The separation of these terms is done for clarity reasons, to distinguish the fixed past values from the decision variables of the optimization problem.

Battery Energy Storage Systems

Finally, the constraints related to the BESS at node $j$ are given as

[TABLE]

where $E_{\textrm{cap,j}}^{\textrm{B}}$ is the installed BESS capacity connected at node $j$; $SoC_{\textrm{min}}^{\textrm{B}}$ and $SoC_{\textrm{max}}^{\textrm{B}}$ are the fixed minimum and maximum per unit limits for the battery state of charge; and, $E_{\textrm{j,t}}^{\textrm{B}}$ is the available energy at node $j$ and time $t$. The initial energy content of the BESS in the first time period is given by E_{\textrm{j,t_{\textrm{0}}}}, and (12c) updates the energy in the storage at each period $t$ based on the BESS efficiency $\eta_{\textrm{B}}$, time interval $\Delta t$ and the charging and discharging power of the BESS $P_{\textrm{j,t}}^{\textrm{B,ch}}$ and $P_{\textrm{j,t}}^{\textrm{B,dis}}$. The charging and discharging powers are defined as positive according to (12d), while (12e) ensures that the BESS is not charging and discharging at the same time. Finally, (12f) limits the reactive power output as a function of the charging or discharging power and the maximum power factor $cos(\phi_{\textrm{max}}^{\textrm{B}})$;

II-B Ancillary services

In grid connected mode, we include the offering of frequency control products following a European market framework [27, 28]. These require power and energy reserves, that can be called at any time. In the following sections, we describe the technical constraints of each product. Please note that only one single FC product is offered at a time, i.e. multiple services are not considered.

II-B1 Primary frequency control (PFC)

PFC is a symmetrical product, i.e. each bid needs to provide symmetrical power bands both for up- and down-regulation, to cover imbalances both from excess production or consumption. The European frequency control reserve cooperation [27] has set the energy requirement to $0.25\cdot P_{\textrm{bid}}$, i.e. the provider has to be able to deliver the full committed power ($P_{\textrm{bid}}$) for a quarter of an hour ($15$ minutes). However, evaluation of realized primary control signals showed that this requirement is conservative[29], i.e. much less energy is actually needed. In this work, only the battery is considered to be able to offer this product. The power reserves for up- and down-regulation are given by

[TABLE]

where $P_{\textrm{bid}}$ is the weekly power size of the PFC bid.

The energy that has to be reserved $\forall t$ is given by

[TABLE]

where $\Delta t_{\textrm{1}}$ is defined to be 15 minutes [30].

II-B2 Secondary Frequency Control (SFC)

SFC is activated after PFC to bring frequency back to the nominal value, and restore the scheduled power exchanges with other control areas. SFC is also symmetrical and requires fast response times. Thus, for the provision of this product, we employ the BESS and the PV units. The power reserves for up- and down-regulation (again symmetrical) are given by

[TABLE]

where $P_{\textrm{bid}}$ is the weekly power size of the SFC bid. Secondary control is activated after a few seconds and is typically completed after 15 minutes [28]. However, in reality this scheme does not guarantee that the energy requirement will not exceed the energy required for a provision of $P_{\textrm{bid}}$ for 15 minutes. By design, there is a continuous secondary call signal that needs to be followed, not accounting for specific energy requirements. For this reason, a statistical approach was followed to analyze ex-post the SFC signal over 1 year in Switzerland, and subsequently derive hourly worst case energy requirements per bid size of secondary frequency power. The worst case values for a 24-h rolling horizon required an energy content of around 5.5 hours times the amount of the bid size in either direction [29]. However, these values are too conservative, and would limit drastically the flexibility on the secondary frequency control market. Thus, we consider as additional constraints only the first $4$ hours of the worst case requirements. Afterwards, the missing/surplus energy can still be bought/sold at the spot market with a lead time of one hour [31].

Furthermore, since PV forecasts are subject to some short-term adjustments, we require that at least 50% of the energy of a worst case call has to come from the BESS.

Thus, the energy content evolution for the first 4 hours of the worst case call is described by

[TABLE]

[TABLE]

where $E_{j,t+\vartheta}^{\textrm{B,2,+}}$ (resp. $E_{j,t+\vartheta}^{\textrm{B,2,-}}$) is the BESS energy content at time t + $\vartheta$ for a call of up (resp. down) regulation at time t; $\vartheta\in\{1,2,3,4\}$ denotes the time for the first 4 hours of the worst case calls, e.g. $E_{j,t}^{\textrm{B,2,+}}$ and $E_{j,t}^{\textrm{B,2,-}}$ correspond to the initial BESS content when the SFC call occurs; and \Delta t_{\textrm{2,\vartheta}}^{\textrm{\pm}} denotes the worst case up- and down-regulation delivery time at hour $\vartheta$ (whereby \Delta t_{\textrm{2,\vartheta}}^{\textrm{\pm}}\leq\Delta t) derived empirically by the ex-post analysis of the SFC signal [29].

The battery energy content for each individual case, i.e. $\forall t,\vartheta$, are required to stay within the acceptable boundaries,

[TABLE]

II-B3 Tertiary Frequency Control (TFC)

Tertiary control is asymmetric (up and down) and significantly slower than PFC and SFC, allowing also flexible loads to participate. For this product, both weekly bids as well as bids for single 4-hour blocks can be provided. In the latter case, the constraints apply only to these 4 hours. The equations are similar to the case of SFC; however, the amount of energy reserves is defined exactly by the regulation of this frequency product, without the need of setting empirical additional constraints. Throughout the duration of the four hours, the full amount of power has to be dispatchable.

Similar to the case of secondary control, a minimum share of energy has to be provided by the BESS. Here, we define that PV generation combined with flexible loads can account for a maximum share of 80% of a call. The power and energy constraints for up-regulation are given by

[TABLE]

[TABLE]

where $P_{\textrm{bid}}$ is the weekly or 4-hour block power size of the TFC bid, and $\Delta t_{\textrm{3}}$ is fixed to 1 hour. The SoC constraint $\forall\vartheta\in\{1,2,3,4\}$ is given by

[TABLE]

The case of down regulation is similar and straightforward. Finally, for all cases the maximum bid size is constrained by

[TABLE]

where $P_{\textrm{t}}^{\textrm{bid,max}}$ is the maximum power size of the FC product. We use the same variable ($P_{\textrm{bid}}$) for the different FC products, because only one can be offered at a time, i.e. we do not consider provision of multiple services by BESS [32].

II-C Islanded mode

In this work, we consider the capability of the distribution grid to be operated in islanded mode, i.e. as a microgrid disconnected from the higher grid level. This is treated by introducing a second set of variables. Most of these constraints are the same as the equations for the grid connected mode and can simply be duplicated.

In this work, the goal in islanded mode is to serve as much of the critical load as possible during the first 24 hours. To achieve that, we utilize the PV generation, BESS and load curtailment. We treat flexible loads as not critical and thus these loads are not considered in the islanded mode. The power balance equations are given by

[TABLE]

where $\alpha_{\textrm{j,t}}^{\textrm{serv,isl}}$ denotes here the fraction of active power served.

Modern grid codes require a minimum power factor requirement in the grid-connected case[26]. However, in the islanded mode we exploit the full functionality of the PV and BESS inverters. Thus, the reactive power provision is described by

[TABLE]

Finally, all the constraints concerning the OLTC are not active in the islanded case. The only link between the set of variables in the grid-connected and the islanded mode is the BESS energy content at timestep $\tau$, when the islanding operation begins, i.e.E_{\textrm{j,\tau}}^{\textrm{B,isl}}=E_{\textrm{j,\tau}}^{\textrm{B}}. After that, the two sets of variables describe independent possible future developments.

III Handling of uncertainty and chance-constrained OPF formulation

This section first describes how the uncertainty is considered in form of chance constraints and then summarizes the final centralized CC-OPF formulation.

III-A Accounting for Uncertainty through Chance Constraints

In order to consider the impact of generation uncertainty, we follow our previous work [23, 33] and we re-formulate the problem using chance constraints [34, 35]. We assume that the PV power injection is the only source of uncertainty (load uncertainty can be also included in a similar way) and we use as input forecast error distributions with different forecasting horizons (1 to 24 hours ahead).

Following [23, 33] we model the voltage and current constraints as chance constraints that will hold with a chosen probability $1-\varepsilon$, where $\varepsilon$ is the acceptable violation probability. E.g., the voltage and current magnitude constraints are reformulated as $\mathbb{P}\left\{V_{\textrm{min}}\leq|V_{\textrm{j,t}}|\leq V_{\textrm{max}}\right\}\geq{1-\varepsilon}$ and $\mathbb{P}\left\{|I_{\textrm{i,t}}^{\textrm{br}}|\leq I_{\textrm{i}^{\textrm{max}}}\right\}\geq{1-\varepsilon}$, respectively. To solve the resulting CC-OPF, we interpret the probabilistic constraints as tightened deterministic versions of the original constraints following the work of [35, 34]. The tightening represents a security margin against uncertainty, i.e., an uncertainty margin. Thus, we express the voltage and current constraints as

[TABLE]

where $\Omega_{\textrm{V}}^{\textrm{lower}},~{}\Omega_{\textrm{V}}^{\textrm{upper}}$ are the tightenings for the lower and upper voltage magnitude constraints and $\Omega_{I_{\textrm{br}}}$ are the tightenings of the current magnitude constraints. The interested reader is referred to [23] for more details on this method.

The uncertainty margins are constant within the OPF solution process, and evaluated outside of the OPF solution. Thus, we use a Monte Carlo approach and the non-linear AC power flow equations to evaluate the boundaries. This further allows us to include any uncertainty probability distribution.

Hence, we form empirical distributions for the voltage and current chance constraints at each time step based on the results from the Monte Carlo simulations. To enforce a chance constraint with $1-\epsilon$ probability we need to ensure that the $1-\epsilon$ quantile of the distribution remains within the bounds. Thus, the tightening corresponds to the difference between the forecasted value with zero forecast error and the $1-\epsilon$ quantile value evaluated based on the empirical distribution resulting from the Monte Carlo Simulations, e.g. $|V_{\textrm{bus,j,t}}^{\textrm{0}}|$ and |V_{\textrm{bus,j,t}}^{\textrm{1-\epsilon\%}}| for the voltage constraints. The empirical uncertainty margins to be used in the next iteration are then given by

[TABLE]

where superscript 0 indicates the current or voltage magnitude at the operating point with zero forecast error. Finally, an iterative algorithm is needed, because the uncertainty margins rely on the derived DER setpoints [36, 34]. Consequently, we alternate between solving a deterministic OPF with tightened constraints, and calculating the uncertainty margins $\Omega_{\textrm{V}}^{\textrm{lower}},~{}\Omega_{\textrm{V}}^{\textrm{upper}},~{}\Omega_{I_{\textrm{br}}}^{\textrm{upper}}$. When the change in the tightening values between two subsequent iterations is below a threshold $(\eta_{V}^{\Omega},~{}\eta_{I}^{\Omega})$, then the algorithm has converged.

III-B Solution Algorithm

In this section, we summarize the proposed solution method for the centralized CC-OPF scheme implemented in an MPC fashion, sketched in Fig. 1. First, the initialization stage sets the uncertainty margins to zero and initializes the voltage levels to a flat voltage profile. At the core of the proposed methodology lies the formulation of the multi-period centralized CC-OPF, which considers the provision of ancillary services as well as the possibility for islanded operation. The CC-OPF calculates the optimal DER setpoints based on a single sweep of the BFS algorithm. The BFS power-flow algorithm then runs until convergence for the obtained control settings. The CC-OPF is then performed again using the updated voltages from the full BFS. These inner iterations are carried out until convergence. After the multi-period OPF has converged, the uncertainty margins are evaluated in the outer loop as described in Section III-A. The iteration index of the OPF loop is denoted by $k$ and the iteration of the uncertainty loop by $m$. The iterative procedure continues until all parts of the algorithm have reached convergence. Then, only the optimal setpoints of the first time step are implemented. Subsequently, the PV forecast is updated, the current timestep is increased and the next CC-OPF problem with a horizon of 24 hours is solved.

The resulting optimization problem is a mixed-integer quadratically constrained program (MIQCP) and can be solved efficiently by modern powerful solvers. The computational burden depends on the dimensions of the grid, the acceptable violation probability, and the number and complexity of the considered DGs.

Due to the efficient handling of the power flow equations through the BFS formulation, hundreds of nodes and branches can be handled without a drastic increase in the computational burden. Regarding the uncertainty handling, the selection of $\epsilon$ influences the execution time of the proposed scheme, since it modifies the feasible area of the optimization problem. The larger the required fulfillment (small values of epsilon), the smaller the feasible area of the optimization problem, making the optimization more demanding. If it is necessary to reduce the computational burden, DGs with complex modes can be handled with reasonable approximations. E.g. constraint (12e) could be replaced as in [22] to avoid the need for binary variables, and the operation of the tap changers could be modeled as continuous variables, rounded ex post to the closest integer.

Overall, however, realistic distribution grid dimensions require solving time in the range of minutes, which is acceptable for such kind of steady state analysis and can be implemented in existing active distribution grids.

IV Case Study - Results

In order to demonstrate the proposed method, we use a typical European radial LV grid [37], sketched in Fig. 2. The installed PV capacity is expressed as a percentage of the total peak load as follows: PV nodes = [12, 16, 18, 19], PV share (%) = [35, 25, 30, 45]. Furthermore, we consider flexible loads up to 5 kW at nodes [17, 18, 19], i.e. $5\%$, $15\%$ and $10\%$ of the corresponding nominal load. The BESS capacity at node 2 is 484 kWh, and the maximum power 484 kW. In this work, we only consider balanced, single-phase system operation, but the framework can be extended to three-phase unbalanced networks as we explain in [24].

The spot market prices were assumed equal to the realized values of 2016 [31]. The realized reserve prices of 2016 are available in [28]. To adjust the cost for the islanded case, we used a constant of $f_{\textrm{isl}}=0.1$, and very high load curtailment cost of $c_{\textrm{t}}^{\textrm{curt,l,isl}}=250\frac{\textrm{€}}{\textrm{MWh}}$. Furthermore, a realized primary control signal was derived from a frequency signal with a temporal resolution of one second. A realized secondary control signal with the same time resolution was taken from [38].

Regarding the uncertainty modeling, we use historical forecast error distributions from an area in Switzerland provided by [38] and we enforce the chance constraints with an $\epsilon=5\%$ violation probability. We assume a maximum acceptable voltage of $1.1$ p.u and cable current magnitude of $1$ p.u. on the cable base. The minimum acceptable voltage is set to $0.9$ p.u..

Using this system, we investigate the capability of the DN to switch to islanded mode, while offering frequency control products. Furthermore, we show how the DN responds to a frequency control call, respecting the islanding requirement. The implementation was done in MATLAB. For the centralized OPF-based control, YALMIP [39] was used as the modeling layer and Gurobi [40] as the solver. The results were obtained on an Intel Core i7-2600 CPU and 16 GB of RAM.

IV-A Islanded operation

The first part of the results refers to the ability of the DN to switch to the islanded mode, where at least 10$\%$ of the load should be served for the next 24 hours. This parameter is estimated to cover emergency services.

IV-A1 Determination of minimal BESS size

A minimum battery energy capacity is required in order to ensure islanded feasibility under different PV injection and loading conditions. Thus, we used historical values of available PV and load data to determine the minimum BESS requirement for islanded operation. We performed yearly MPC-OPF calculations with a 24-hour horizon, without considering uncertainties, to estimate the needed BESS size iteratively; i.e. we kept increasing the BESS size until we derived feasible solutions for the whole year. The worst case period is shown in Figure 3, indicating a minimum BESS size of $220$ kWh.

In order to allow provision of AS, we investigated various BESS capacities corresponding to $1.4-2.6$ times the needed minimum value. In the remaining simulations, we will consider a BESS of $484$ kWh.

IV-A2 Switch to islanded mode

According to Section IV-A, the switch to islanded operation should be feasible at any time instant. Figure 4 shows the evolution of the BESS SOC for islanding at distinct hours in the considered time period. The power balance is kept using the BESS capacity, PV injections, and load and PV curtailment. As can be observed, the BESS SoC evolution depends on the PV generation and load forecasts; At noon hours, the PV units provide power for the loads and BESS charging, while at night the BESS is discharged to guarantee a 24-hour islanded operation.

IV-B Frequency control

In the grid-connected case, the DN offers frequency regulation as an ancillary service, while at the same time fulfilling the islanded requirement for the next 24 hours.

IV-B1 PFC

Assuming that the BESS is always charging or discharging at a maximum rate of 1C111A C-rate is a measure of the rate at which a BESS is charged or discharged relative to its maximum capacity. A 1C rate means that the discharge current will discharge the entire battery in 1 hour. to limit the capacity fading from offering frequency control products [41], an energy requirement of 15 minutes PFC power in both directions, i.e. 30 minutes in total, translates into reserving 50% of the total BESS storage capacity. The algorithm keeps the SoC at the upper limit to minimize load shedding in case of a switch to the islanded mode.

Figure 5 shows the BESS SoC while providing PFC over a summer week. Staying outside of the red area guarantees that in the case of a switch to islanded operation at any time step the critical load can be supplied by preserving a minimum BESS energy content based on load and PV generation forecasts. The orange area represents the energy limit imposed by the offered frequency control product. The white area defines the allowable feasible region for the SOC, with the black line showing the optimization result. In case of overlapping between the orange and red area, the more limiting area is relevant. In case of operating in islanded mode, frequency reserves are not provided anymore.

IV-B2 SFC

For this product we consider also PV units, which can curtail power providing down-regulation. Hence, the upper bound on the energy level of the storage during hours with PV injections is relaxed, as seen in Fig. 6. The BESS can be charged during these hours, leading to higher self-consumption and more available energy in case of a switch to islanded mode.

IV-B3 TFC - weekly offer

Up regulation

Providing maximum up TFC regulation resulted in a fully charged BESS, as can be observed in Fig. 7. In this way, we not only achieve maximum reserve provision, but also minimum load curtailment in the islanded mode. Limited flexibility is offered by flexible loads, as can be seen by the white areas, the size of which does not influence the maximum bid size.

Down regulation

The case of down regulation is shown in Fig. 8, where the optimization tries to keep the SoC low in order to respond to a TFC dispatch call, while at the same time respecting the islanding requirement. Similar to the SFC case, during noon hours with solar power, the SoC can be increased, since PV power curtailment is available.

IV-C Call for SFC

So far, we studied the needed power and energy reserves. In this section, we simulate the response of the DN to an actual continuous SFC signal. Since we cannot forecast the signal, we used the realized signal from 2016. Figure 9 shows the worst-case week in terms of needed power of the SFC signal as well as the corresponding cumulative energy requirement. We consider the possibility of participating in the spot market with a lead time of four hours. As can be observed, the algorithm chooses to buy energy on the spot market three times indicated by the red circles. The dashed line corresponds to the cumulative energy demand without spot market purchases, whereas the solid line to the resulting cumulative energy given the purchases in the spot market. Finally, Fig. 10 shows the evolution of the BESS SoC following the SFC signal in solid, and the SoC without offering SFC with a dashed line. As can be observed, the three purchases of power are needed so that the SoC is kept high enough to allow for the islanded mode.

IV-D Impact of BESS size on the rating of the transformer

As a final case study, we investigate the impact of the BESS size on the needed rating of the MV/LV transformer, without offering frequency control products. The BESS can contribute to the power needed to and from the active distribution grid, reducing the required transformer rating. In this way, the service of investment deferral can be offered to the operator, which might need to cope with increasing demand or DG injections.

Figure 11 shows the required rating of the secondary substation transformer, varying the energy capacity of the BESS placed at the same node. A seasonal analysis allows calculating the most critical period, i.e. winter in our case, that defines the needed rating. We observe that the larger the energy BESS capacity, the smaller the required transformer rating; however, the BESS contribution is decreasing with increasing BESS size.

V Conclusion

Modern DNs consider the active control capabilities of DERs in order to provide a secure, reliable and optimal operation of the grid. Furthermore, they can offer ancillary services to higher voltage levels, or even operate disconnected from the main grid.

In this paper, we have shown that ADGs can be coordinated through centralized control schemes to provide ancillary services and provision for islanded operation. The proposed method allows ADGs to support the transmission network but at the same time provide increased resilience through controlled islanding. We have shown how the different operational requirements can be formulated in the problem constraints and provided techniques to tackle the uncertainty.

Acknowledgments

The work of P. Aristidou was partially supported by the Engineering and Physical Sciences Research Council (EPSRC) in the UK under grant reference EP/R030243/1.