Provision of Primary Frequency Response as Ancillary Service From Active Distribution Networks to the Transmission System

IEEE TRANSACTIONS ON SMART GRID, VOL. 12, NO. 6, NOVEMBER 2021 4971 Provision of Primary Frequency Response as Ancillary Service From Active Distribution Networks to the Transmission System Eleftherios O. Kontis ,Member, IEEE, Álvaro Rodríguez del Nozal , Juan M. Mauricio, Senior Member, IEEE, and Charis S. Demoulias ,Senior Member, IEEE Abstract—This paper deals with the provision of primary frequency response (PFR) as ancillary service (AS) from active distribution networks (ADNs) to the transmission system (TS). In particular, two methodologies are developed. The first one aims to quantify the PFR capability range of the ADN. This range is defined by determining the range of the aggregated, i.e., equivalent, active power - frequency P(f)droop curves that can be provided at the point of interconnection (POI) with the TS. The second one targets to optimally control P(f)droop curves of individual distributed energy resources (DERs), installed in the premises of the ADN, to guarantee specific frequency regulation characteristic at the POI. This frequency regulation characteristic is expressed by means of a P(f)droop curve. Both methods are tested on two discrete distribution systems. Several test cases are examined to demonstrate their implementation. Additionally, comparisons against conventional approaches and time series simulations are conducted to evaluate the performance of the proposed methods. Index Terms—Active distribution networks, ancillary services, distributed energy resources, primary frequency response. I. INTRODUCTION TRADITIONALLY, after a large generation or load disturbance, power system frequency is controlled by utilizing inertial and primary frequency response (PFR) of synchronous generators (SGs) [1]. However, nowadays, power system industry is moving towards generating electricity from distributed energy resources (DERs), that are mainly connected to the utility grid via power electronic converters [2]. Consequently, conventional SG-based power systems are transformed into converter-based ones [3]. This shift Manuscript received January 8, 2021; revised April 23, 2021 and July 19, 2021; accepted July 28, 2021. Date of publication August 6, 2021; date of current version October 21, 2021. This work was supported by the European Union, Horizon 2020 Project “EASY-RES” under Grant 764090. Paper no. TSG-00044-2021. (Corresponding author: Eleftherios O. Kontis.) Eleftherios O. Kontis is with the School of Electrical and Computer Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece, and also with the Department of Electrical and Computer Engineering, University of Western Macedonia, 50100 Kozani, Greece (e-mail: ekon- [email protected]). Álvaro Rodríguez del Nozal and Juan M. Mauricio are with the Department of Electrical Engineering, Universidad de Sevilla, 41092 Seville, Spain. Charis S. Demoulias is with the School of Electrical and Computer Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece. Color versions of one or more figures in this article are available at https://doi.org/10.1109/TSG.2021.3103060. Digital Object Identifier 10.1109/TSG.2021.3103060 towards converter-interfaced generation modifies the dynamic properties of power systems [4], affecting also their frequency control capability and performance [5]. Indeed, converter-interfaced units do not possess inherent inertia and primary frequency control capabilities [3], [5]. Thus, significant frequency changes, fast-frequency dynamics, under frequency load shedding and even frequency instabilities are expected in power systems dominated by converter-interfaced units [5], [6]. Therefore, to enhance system stability and to ensure the uninterruptible provision of high quality power, frequency control functionalities, such as inertia and PFR, should be embedded to converter-interfaced DERs [4], [5]. In fact, nowadays several interconnection standards impose that DERs should provide frequency support to the power system. Methodologies to provide frequency control functionalities from DERs are reviewed in [7] and [8], while several standards discussing frequency control requirements for the interconnection of DERs are summarized in [4]. The above-mentioned frequency control functionalities allow large-scale converter-interfaced power plants, that are directly connected to the transmission system (TS), to participate in balancing and ancillary services (ASs) markets [9]. However, nowadays, a large amount of green energy is also produced by small-scale DERs, that are connected to active distribution networks (ADNs) [10]. Therefore, as discussed in [11], equal opportunities should also be provided to these units to increase their value and revenues by participating in the above-mentioned AS markets. In this context, several methods have been proposed during the last years to control the operation of DERs, connected to ADNs, in order to provide AS to the TS. For instance, in [12]–[14] methods to quantify the active and reactive flexibility area of ADNs are developed. Congestion management of TSs through the utilization of ADN assets is suggested in [15], while in [16] coordination of DERs is proposed to provide voltage support to TSs. Given the focus of this work, we provide in the next paragraphs a review of existing methods aiming to coordinate the operation of several DERs in order to provide PFR as AS to the TS. In [17]–[19] centralized control architectures are developed to facilitate the provision of PFR from large-scale wind farms (WFs) to the TS. In particular, in these papers, optimization problems are formulated aiming to optimally distribute PFR among the wind turbines of a WF based on different This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License. For more information, see https://creativecommons.org/licenses/by-nc-nd/4.0/ 4972 IEEE TRANSACTIONS ON SMART GRID, VOL. 12, NO. 6, NOVEMBER 2021 operational objectives, e.g., cost minimization [17], power loss minimization [18], maximization of kinetic energy [19], etc. The concept of virtual power plant (VPP) is introduced in [20] to coordinate the operation of combined wind and photovoltaic (PV) power plants for the provision of frequency ASs to the TS. In [21] a VPP consisting of WFs and demand response (DR) units is proposed, targeting to aggregate the contribution of several resources in order to provide PFR as AS to the TS. The aggregator concept is utilized on [22] to optimally coordinate the operation of multiple battery energy storage systems (BESSs). However, all the above-mentioned methods are oriented to large-scale DERs and BESSs, that are directly connected to the TS. As a consequence, they cannot be applied for ADN analysis, since they neglect load peculiarities and ADN characteristics. To overcome these issues, methods oriented on ADN and microgrid (MG) analysis are developed in [3], [23]–[26]. In these methods, the aggregator concept is utilized to coordinate the operation of several ADN assets, such as electric-vehicles (EVs), PVs, DR units, and BESSs, in order to provide in an aggregated manner frequency ASs to the utility grid. The coordination of ADN assets is performed via optimization routines. However, neither of the above-mentioned works integrates the network modeling within the optimization. Thus, technical and operational constraints of ADNs, such as grid voltage and line thermal limits, are not fully considered [27]. In [28] a stochastic frequency control method for gridconnected MGs, operated under active power - frequency P(f) droop curves, is proposed. However, MGs are simulated using aggregated reduced order models. Thus, local issues such as voltage violations and congestion problems are not considered. The provision of PFR from EVs is discussed in [29] and [30]. In particular, in these works algorithms to determine charging/discharging power of EV batteries are formulated, aiming to control the power delivered at the point of interconnection (POI) with the TS. However, in these approaches line thermal limits are not considered. Another approach is presented in [31] and [32]. Specifically, in these works, P(f)droop curves of individual DERs, installed within an ADN, are designed to achieve a specific PFR at the POI with the TS. To determine the individual P(f)droops, the contribution of each DER to the aggregated PFR is evaluated using participation factors, derived via linearized power flows. However, these methods are designed based on the assumption that all DERs contribute to the PFR proportionally to their size. Having determine the gaps in the existing literature, scope of the paper is to develop tools and methods in order to facilitate the provision of PFR from ADNs to the TS. For this purpose, a methodology is initially developed to accurately determine the PFR capability range of ADNs, i.e., the minimum and maximum PFR that can be offered at the POI with the TS. To achieve this, the aggregated contribution of all DERs to PFR is quantified by determining the aggregated, i.e., equivalent, P(f)droop curves that can be provided at the POI with the TS. Additionally, a second methodology is proposed to optimally coordinate the P(f)droop curves of individual DERs in order to provide specific PFR characteristics at the POI with the TS. The dispatch of DERs is performed in an optimal way, targeting to minimize an objective function. Both methodologies are developed based on conventional optimal power flow (OPF) formulations, incorporating this way ADN modeling as well as grid technical and operational constraints into the optimization problem. Thus, contrary to the existing approaches, the proposed methodologies fully exploit operational limits of ADNs, maximising their contribution to PFR. Additionally, the use of conventional OPF formulations facilitates the integration of the proposed methods to commercial power system software analysis tools. To the best of the authors’ knowledge, no similar methodologies are reported in the literature. Finally, a framework is developed to define work flows, signal and information exchange between TS operators (TSOs) and distribution system operators (DSOs), that are responsible for the operation of ADNs. The rest of the paper is organized as follows: Section II presents the theoretical background, explaining the operation of DERs under the P(f)droop control scheme and discussing the technical and operational constraints of ADNs. The main idea of providing PFR from ADNs to the TS is explained in Section III. Section IV presents the methodology developed to estimate the PFR capability range of ADNs, while in Section V the method proposed to optimally dispatch DERs is presented. In Section VI, the proposed methodologies are tested on a single feeder distribution system. Comparisons against conventional approaches are carried out and the impact of forecast errors on the accuracy of the proposed methodologies is quantified. More evaluation results are provided in Section VII, where both methodologies are tested on a double feeder distribution system. In Section VIII a discussion concerning the applicability of the proposed methods is performed. Additionally, technical and regulatory issues that must be addressed to allow the implementation of the proposed methods are discussed. Finally, Section IX summarizes the research findings and concludes the paper. II. THEORETICAL BACKGROUND Operation of DERs under P(f)droop control and modeling of ADNs are briefly reviewed in this section to provide the required background for the analysis of the proposed methodologies. A. Operation of DERs Under the P(f) Droop Control To participate in PFR, DERs should be deloaded, i.e., operated with a headroom under their maximum point power (MPP), to ensure that sufficient generation margin is available at any time instant [33]. Methodologies to design de-loading controllers for wind turbines and PV systems are summarized and reviewed in [7] and [8]. Additionally, DERs should be controlled with P(f)droop curves to adjust their active power based on grid frequency. A typical P(f)droop is presented in Fig. 1. Based on grid frequency f, three regions (modes of operation) are identified. Normal operation mode:fn−f 2≤f≤fn+f 2. Grid frequency is around its nominal value, i.e., fn, within a predetermined deadband f.fis determined by the TSO. DERs inject an amount of active power which is frequency insensitive and is KONTIS et al.: PROVISION OF PFR AS AS FROM ADNS TO TS 4973 Fig. 1. Typical P(f)droop curve. defined as Pk≡Pc=PMPP,k(1−ph,k). (1) PMPP,kis the maximum available power of the k-th DER; ph,kis expressed in p.u. and denotes the headroom, i.e., the available power reserve, of the k-th DER. Under-frequency mode:fmin ≤f<fn−f 2.Herefmin is the minimum permissible grid frequency, determined by the TSO. Active power injection of DERs is defined as Pk=PMPP,k+PMPP,kph,k(f−fmin) fmin +f 2−fn .(2) Over-frequency mode:fn+f 2<f≤fmax.Herefmax is the maximum permissible grid frequency, determined by the TSO. Active power injection of DERs is defined as Pk= Pc(f−fmax)+Pmin,kfn+f 2−f fn+f 2−fmax .(3) Pmin,kis the minimum power that k-th DER can inject. B. Technical and Operational Constraints of ADNs ADNs are connected to the TS via interconnection transformers (ITFs). For the modeling of ADNs, in this paper, an undirected graph G=(V,E)is used. V={0,1,...,p}is the set of ADN buses. Without loss of generality, it is assumed that the ITF is placed between buses 0 and 1. Bus 0 denotes the high voltage (HV) side of the ITF, while bus 1 the MV side. The set of power lines is represented using E⊆V×V. Moreover, Ni:={j:(j,i), (i,j)∈E}is used to represent the neighborhood of bus i, i.e., the set of buses that are connected to bus i. Finally, Di={1,2,...,di}is used to denote the set of DERs that are installed at bus i. To simulate ADNs, the following set of equations is used  k∈Di Sk−SL i=√3UiI∗ i,∀i∈V(4)  j∈Ni Ii,j=Ii,∀i∈V,j∈Ni(5) Ui−Uj=Zi,jIi,j,∀i∈V,j∈Ni.(6) SL iis the complex load consumption at bus i.Skis the complex power of the k-th DER at bus i. Real power injection of the k-th DER depends on grid frequency and computed using (1) -(3).Uiand Iidenote the complex phase-to-phase voltage and the complex current injection at bus i, respectively. Ii,jis the current flowing through line (i,j);Zi,jis the corresponding line impedance. Note that Ii,j=−Ij,iand Zi,j=Zj,i. To guarantee the secure operation of the ADN, line currents and bus voltages should remain within the following limits:  Ii,j ≤Imax i,j,∀i∈V,j∈Ni(7) Vmin ≤|Ui|≤Vmax,∀i∈V(8) here Imax i,jis the maximum ampacity of the line, while Vmin and Vmax are the lower and upper voltage limits, respectively, specified by relevant grid codes and power quality standards such as EN50160. III. PROVISION OF PFR FROM ADNSTOTHETS To provide PFR as AS to the TS, DSOs shall coordinate the operation of several DERs. In particular, DSOs should perform the following tasks. i) Estimate the minimum and maximum PFR capability of ADNs. Specifically, DSOs should evaluate the range of the aggregated, i.e., equivalent, P(f)droop curves that can be provided at the POI with the TS. To evaluate this range, technical constraints must be considered to ensure the safe grid operation. Information concerning PFR range will allow DSOs to actively participate in future AS markets [32], [34]. ii) Provide specific, TSO-defined, PFR. In future AS markets, ADNs will compete with other market participants concerning the provision of PFR as an AS to the TS. All market offers will be evaluated by the market operator targeting to minimise the cost of the AS. Thus, during specific time instants ADNs will have to reduce their market share, simply because other market participants will provide more competitive offers/products. In these cases, TSOs should send appropriate dispatch signals to the corresponding ADNs to reduce their contribution to PFR, since they will not be remunerated for the extra power. Towards this objective, TSOs should determine PFR characteristics of individual ADNs by defining the P(f)droop curve that each ADN must provide at the POI with the TS. The DSO should then dispatch DERs of the ADN to achieve this TSO-defined P(f)droop, ensuring also that all ADN technical constraints are fulfilled. Dispatch must be performed in an optimal way to minimize the cost of the AS. To implement the above-mentioned tasks, the framework of Fig. 2 is proposed. This framework describes the required work flows as well as signal and information exchange between TSOs and DSOs in order to facilitate the provision of PFR from ADNs to the TS. As shown, during predetermined time intervals (t), TSO informs the DSO concerning permissible frequency limits, i.e., f,fmax, and fmin. Generally, t is a TSO-defined parameter that can vary from a few minutes up to several hours. Nevertheless, according to [35], each 4974 IEEE TRANSACTIONS ON SMART GRID, VOL. 12, NO. 6, NOVEMBER 2021 Fig. 2. Interaction between DSO and TSO for the provision of PFR. Summary of paper contributions. TSO should perform contingency analysis at least once every 15 minutes. Thus, in this paper, tis set to 15 minutes. Additionally, f,fmax, and fmin limits are generally defined by interconnection standards [4]. However, their values can be modified, practically on a monthly basis, by TSOs [35]. Thus, a regular update is required. In this paper, this update is performed every t. Using f,fmax, and fmin, DSO evaluates the PFR capability range of the ADN. i.e., the range of the equivalent P(f) droop curves that can be provided at the POI with the TS. For this purpose, the maximum (Pmax 0) and minimum (Pmin 0) power that the ADN can deliver at the POI are computed. Note that Pmax 0and Pmin 0are evaluated for several grid frequencies varying from fmin to fmax. This approach is proposed since load consumption is affected by the grid frequency [36]. PFR capability range is eventually defined by Pmax 0,Pmin 0,fmax and fmin limits. An example is given in Fig. 3. In this case, PFR range is determined by the hatched area. Then, the PFR capability range is forwarded to the TSO, which conducts a contingency analysis to determine the overall required PFR reserves. Contingency analysis is performed based on normative incidents. An example of how contingency analysis should be conducted is presented in [37]. Afterwards, the contribution of each ADN to PFR is determined through market mechanisms. Subsequently, based on the contribution of each ADN, the TSO defines the corresponding P(f)droop curve and sends appropriate set points to the DSO. These set points determine the power that the ADN must provide at the POI with the TS when grid frequency is equal to fmax,fn, and fmin. For instance, to define the P(f)droop that is plotted with Fig. 3. Range of aggregated P(f)droop curves at POI. Line connecting points Aand Ddenotes Pmin 0limits. Line connecting points Band Cdenotes Pmax 0limits. The blue curve denotes a TSO-defined P(f)droop curve. the blue curve in Fig. 3, the TSO should send set points E (P1,fmax), F(P2,fn), and G(P3,fmin) to the DSO. In the final stage of the proposed framework, DSO optimally dispatches DERs to acquire the TSO-defined P(f)droop curve at the POI. The dispatch is performed through an optimization routine, targeting to minimize the cost of the AS. The algorithm proposed to estimate PFR capability range is presented in Section IV. Moreover, the algorithm proposed to optimally dispatch DERs to guarantee a specific TSO-defined P(f)droop curve at the POI is presented in Section V. IV. EVALUATION OF THE PFR CAPABILITY To estimate the PFR capability range of ADNs, the algorithm of Fig. 4a is developed, consisting of seven main steps. Step 1 (Initialization phase): The topology of the ADN, the location of DERs and loads as well as short-term generation and load forecasts are provided as inputs to the proposed algorithm. Additionally, DSO specifies grid technical constraints. These constraints refer to current and voltage limits and they are expressed using (7) and (8), respectively. Step 2 (Derivation of ADN model): The ADN model is developed for grid frequency f. As shown, in the first iteration of the algorithm, all ADN parameters are calculated for fmin. Step 3 (Determination of the maximum power that can be delivered at POI): An OPF is performed. Aim of the OPF is to adjust power injections of all DERs, i.e., Pkand Qk, in order to compute the maximum power Pmax 0that can be provided under grid frequency fat the POI. The OPF is expressed by (4) - (8). The objective function of the OPF and the associated control variables are formally defined as ¯ Pmax,¯ Qmax=arg max ¯ P,¯ QP0.(9) Here, ¯ P={Pk,∀k∈Di,∀i∈V}and ¯ Q={Qk,∀k∈Di,∀i∈ V}declare the sets of active and reactive powers of all DERs, respectively. Super-index max denotes the value of those power injections that maximize (9). P0is the active power flowing through the POI. KONTIS et al.: PROVISION OF PFR AS AS FROM ADNS TO TS 4975 Fig. 4. a) Algorithm-1: Methodology to estimate the PFR capability of ADNs. b) Algorithm-2: Methodology to dispatch DERs in order to provide specific TSO-defined PFR at POI. Step 4 (Determination of the minimum power that can be delivered at POI): A second OPF is conducted. The objective of this OPF is to adjust power injections of all DERs in order to determine the minimum power Pmin 0at POI. As previously, OPF is expressed by (4) - (8). However, in this case, the objective function is ¯ Pmin,¯ Qmin=arg min ¯ P,¯ QP0.(10) Super-index min denotes power injections that minimize (10). Step 5 (Save results and proceed to the next frequency): At this step, grid frequency fand the corresponding Pmax 0and Pmin 0values are stored in three separate vectors. Afterwards, grid frequency is increased by a predetermined step, i.e., df. Step 6 (Condition for algorithm termination): If the resulting frequency is higher than fmax, i.e., the maximum permissible grid frequency, the algorithm proceeds to step 7. Otherwise, as shown in Fig. 4a, the algorithm moves back to Step 2 and the ADN model for the new grid frequency is developed. Step 7 (Forward potential droop range to the TSO): In this step, the algorithm terminates and the vectors containing Pmax 0and Pmin 0values for all considered grid frequencies are forwarded to the TSO. To facilitate the reading, vectors containing Pmax 0and Pmin 0values are depicted in Fig. 3 using purple and orange lines, respectively. These vectors alongside with fmax and fmin limits (red and yellow dotted lines of Fig. 3, respectively) define PFR capability range of the ADN. V. OPTIMAL PFR DISTRIBUTION As already discussed, TSO determines, based on contingency analysis and market mechanisms, the P(f)droop curve that the ADN shall provide at the POI with the TS. TSO communicates the required P(f)droop curve to the DSO by sending appropriate set points. These set points include the power that the ADN should provide at the POI under fmin,fn, and fmax. For the rest of the paper, these powers are denoted as Pfmin,Pfn, and Pfmax , respectively. Subsequently, DSO coordinates the P(f)droop curves of individual DERs in order to guarantee this TSO-defined droop. To coordinate the droop curves of individual DERs, the following solution is proposed: Power injections of each DER, i.e., Pkand Qk, are computed using OPF calculations. Note that the power delivered at the POI changes based on grid frequency, i.e., based on the TSO-defined P(f)droop curve. Thus, Pkand Qkmust be determined for different grid frequencies to ensure that in each frequency the aggregated contribution of all DERs provides the TSO-defined P(f)droop curve. In the simplest case, Pkand Qkshould be determined at least for fmin,fn, and fmax, i.e., for grid frequencies that belong to different regions of the TSO-defined P(f)droop curve. Algorithmic details are presented in Fig. 4b. A stepby-step analysis of the proposed solution is provided below. Step 1 (Initialization phase): Grid topology, location of DERs and loads as well as short-term generation and load foracasts are provided as inputs to the algorithm. Additionally, DSO specifies ADN current and voltage limits using (7) and (8), respectively. Finally, DSO determines the values of weighting factors wloss line and wd k. Further information for these weighting factors are provided in the next Step. Step 2 (Network optimization): Pkand Qkof all DERs are computed for different grid frequencies, varying from fmin up to fmax. At each frequency, Pkand Qkare computed by solving an OPF, targeting to minimize the cost of the AS. The proposed OPF is expressed by (4) - (8). The objective function and the control variables are formally defined as ¯ PPFR,¯ QPFR=arg min ¯ P,¯ QC, s.t. (11) P0=ct,(12) with C=Cline loss +Cd.(13) Here, super-index PFR denotes the values of those power injections that minimize (11). P0is the active power that must be provided at the POI under a specific grid frequency. Cline loss and Cddenote the total cost related with line losses and the remuneration of DERs, respectively. A detailed description is provided below. Line losses: Lines losses are computed as Pline loss = (i,j)∈E Ri,jIi,j2.(14) Here, Ri,jis the resistance of line (i,j). The total cost associated with line losses is evaluated in this paper using the following equation Cline loss =wline lossPline loss,(15) 4976 IEEE TRANSACTIONS ON SMART GRID, VOL. 12, NO. 6, NOVEMBER 2021 where wline loss is a weighting coefficient that penalizes the power losses at the ADN. This weighting coefficient is determined by the DSO. Remuneration of DERs: To provide PFR, DERs should be de-loaded, i.e., operated under their MPP. Inevitably, the deloaded operation of DERs results in loss of green power, that subsequently can lead to loss of revenues for DER owners. To avoid this loss of revenues, DER owners should be remunerated appropriately for the provision of the corresponding AS. The cost function adopted in this paper to define remuneration costs is the following: Cd= k∈Di,i∈V wd kPMPP,k−Pk.(16) Here, wd kis a weighing coefficient, determined by the DSO, that is used to penalize de-loading of DERs. As already discussed, to ensure that the TSO-defined P(f) droop curve is achieved with high accuracy, Pkand Qkmust be computed for several grid frequencies varying between fmin and fmax. As shown in Fig. 4b, this is achieved by iteratively solving (11) - (13) for different grid frequencies. The maximum number of iterations is denoted as m. However, as the number of iterations increases, the total problem complexity is increased as well. Thus, to preserve low computational complexity, injections Pkand Qkare determined in this paper only for three grid frequencies (m=3), belonging to different regions of the P(f)droop. As demonstrated in Section VI-B, this simplification leads to very accurate results. Specifically, in the first iteration, i.e., m=1, grid frequency is equal to fmin and DERs are operated in under-frequency mode. In this case, power at POI, i.e., P0, is constant and equal to Pfmin .In the second iteration, grid frequency is set to fn. DERs operate in normal mode and P0is set to Pfn. Finally, in the last iteration, grid frequency is equal to fmax and DERs operate in over-frequency mode. In this case, P0is equal to Pfmax .As shown in Fig. 4b, for each grid frequency, the corresponding ADN model is developed and a dedicated OPF is solved to determine Pkand Qk. Step 3 (Condition for algorithm termination): If mis higher than m, the algorithm moves to Step 4. Otherwise, the algorithm goes back to Step 2 and (11) - (13) are evaluated for the new grid frequency. Step 4 (Design of individual P(f) droop curves): When all iterations are finished, individual P(f)droop curves are designed, using (17) - (19), and sent to the DERs. Pk≡Pfn,k=Pfmax,k1−ph,k(17) Pk=Pfmax,k+Pfmax,kph,k(f−fmin) fmin +f 2−fn (18) Pk= Pfn,k(f−fmax)+Pfmin,kfn+f 2−f fn+f 2−fmax (19) Here, Pfn,k,Pfmax,k, and Pfmin,kdenote the active power that the k-th DER inject at fn,fmax, and fmin, respectively. ph,k is computed via (17). Note that contrary to the conventional implementation, presented in (1) - (3), in the proposed Fig. 5. System under study. approach DERs do not inject their MPP at fmin. On the contrary they inject Pfmax,k, which can be equal or less to the corresponding MPP. Additionally, DERs do not inject their minimum power at fmax. They inject Pfmin,k, i.e., more or equal power to the corresponding Pmin. This approach is proposed to offer higher flexibility to the TSO and further justified in the next paragraph. TSO may request from the ADN to provide at fmin apower at POI lower than the maximum available. Similarly, TSO may request a power at fmax higher than the minimum available. For instance, TSO may request the blue P(f)droop curve of Fig. 3. To offer this flexibility to TSO, DERs should operate under their MPP at fmin and above their minimum permissible power at fmax. VI. SIMULATION RESULTS ON A SINGLE FEEDER DISTRIBUTION SYSTEM In this section, the proposed methodologies are evaluated on a single feeder distribution system. In particular, the analysis is conducted using the European MV benchmark grid of CIGRE [38]. The required OPFs are solved in Python using Sequential Least Squares Programming (SLSQP) through the minimize function of Scipy library [39]. A. System Under Study The topology of the examined test system as well as the rated power (Pr,k) of DERs and loads are presented in Fig. 5. The nominal grid voltage and frequency are 20 kV and 50 Hz, respectively. For all loads a power factor equal to 0.9 is considered. All loads are modeled as constant power loads. Minimum permissible power of all DERs is assumed equal to 10% of the Pr,kvalue. fmax and fmin are assumed equal to 52 Hz and 46 Hz, respectively, and fis set to 0.4 Hz. Maximum ampacity of power lines is 416 A. Vmin and Vmax are considered equal to 0.95 p.u. and 1.05 p.u., respectively. All other grid parameters are from [40]. B. Explanatory Examples In this section, the functionalities of the proposed methods are demonstrated via indicative examples. To create an easily KONTIS et al.: PROVISION OF PFR AS AS FROM ADNS TO TS 4977 Fig. 6. TSO-defined P(f)droop curve at POI. The power at POI is equal to 12 MW, 10 MW and 8 MW under fmin,fn,andfmax, respectively. The approximation error for different values of mis presented on the right y-axis. observable test case, only the WTs of Fig. 5 are considered in the analysis. All WTs inject their rated power. No loads are considered in the grid. Initially, the PFR capability range of the ADN is evaluated using Algorithm-1. In particular, Pmax 0and Pmin 0 are computed for both fmin and fmax. The resulted PFR range is eventually defined by an area similar to the one presented in Fig. 3, with points A(1.694 MW, 52 Hz), B(14.398 MW, 52 Hz), C(14.400 MW, 46 Hz), and D(1.694 MW, 46 Hz). The analysis reveals that Pmax 0is almost 15% lower compared to the total power of the WTs. Indeed, the proposed method restricts the contribution of WTs to ensure that line thermal limits are not violated. Note that a further increase of Pmax 0will lead to congestion in line connecting nodes N1 and N2. According to the framework of Fig. 2, once the PFR capability range is evaluated, TSO conducts a contingency analysis to determine the P(f)droop curve that the ADN should provide at the POI. However, as already discussed, in this paper we do not propose a new method for contingency analysis. Thus, it is assumed that the TSO conducts the required contingency analysis using his own tools [37]. In order to demonstrate Algorithm-2 we suppose that the TSO requests at the POI the P(f)droop curve of Fig. 6. This TSO-defined droop curve can be obtained using several dispatch scenarios. To display the generic nature of the proposed methodology, two indicative scenarios are presented below. •S1: All WTs participate in PFR. Cost related to line power losses is considered in the optimization problem. To simulate this scenario, the following weighting factors are considered: wline loss =1 and wd k=0 for every k∈Di,i∈V. In this scenario, the dispatch is practically performed in order to minimize the total power system losses. •S2: WT4 does not participate in PFR. The cost of line power losses is not considered in the optimization problem. The weighting factors are: wline loss =0, wd GWT4=1, and wd k=0 for every k∈{GWT2,GWT3,GWT5}. In this scenario, the dispatch is performed targeting to minimize the total power that de-loaded from the WTs, ensuring also that WT4 does not contribute in the PFR. The active power injections of each DER for each scenario are presented in Table I. Power that is de-loaded from each DER as well as system total power losses are shown in Table II. In both scenarios, the power delivered at the POI for fmin,fnand fmax is the one requested by the TSO, i.e., the black TABLE I RESULTS OF THE PFR DISTRIBUTION PROBLEM TABLE II TOTAL SYSTEM LOSSES AND POWER THAT ISDE-LOADED curve of Fig. 6. In S1, line losses are minimized. Consequently the total power provided by the DERs is lower compared to S2. In S2, the main objective is to operate DER GWT4pro- viding as much power as possible. Therefore, GWT4 does not participate in PFR, injecting always its MPP. Power injections of the remaining DERs are adjusted accordingly to provide the desired droop curve at the POI. The impact of the number of iterations (m) on the accuracy of Algorithm-2 is also evaluated, using the assumptions of S2. Results are depicted in Fig. 6. As shown, as the number of iteration increases the absolute error decreases. The maximum error is observed for m=3 and it is equal to 3 kW. This error is merely 0.025% of the power flowing through the POI. Nevertheless, as the number of iterations increases, the required execution time is increased as well. In particular, for mequal to 3, 5, 7, 9 and 11 the execution time is equal to 4.47 s, 8.50 s, 11.20 s, 14.12 s and 17.27 s, respectively. Simulations have been conducted using an i5-7200, 2.71 GHz, 8 Gb RAM PC. For the rest of the analysis, missetto3,since this value achieves a good compromise between computational burden and active power error at POI. C. Comparative Analysis In this section, the methodology proposed to dispatch DERs is compared with the approach developed in [31]. In this work the contribution of each DER to PFR is proportional to the rated power of the unit. Additionally, to allocate the PFR among the installed DERs, the concept of participation factors is introduced. In particular, for each DER a participation factor is computed using sensitivity analysis. Participation factors aim to quantify the impact of individual DERs on the active power flowing through the POI. 4978 IEEE TRANSACTIONS ON SMART GRID, VOL. 12, NO. 6, NOVEMBER 2021 Fig. 7. Comparison with a conventional approach. The TSO-defined P(f)droop curve is depicted in Fig. 7 with black color. To achieve this droop curve at the POI, the P(f)droop curve of individual DERs are designed using both methods. To achieve a fair comparison, the contribution of each DER to the PFR is assumed for both methods proportional to the corresponding rated power. The proposed method is executed assuming (m=3). Concerning the second method, the participation factors are computed as proposed in [31]. The resulted error from both methods is presented on the secondary y-axis of Fig. 7. The error of [31] reaches to a maximum value equal to 7.75 kW. This error is almost four times higher compared to the error of the proposed method. The computed errors for both methods are trivial compared to the total power flowing through the POI, validating their accuracy. However, it is worth noticing that the proposed methodology has a more generic nature compared to the method of [31]. Indeed, in [31] the PFR is proportionally distributed to the installed DERs based on their rated power. On the contrary, using the proposed methodology different scenarios can be considered, e.g., some DERs do not participate at all in the PFR. This is also verified by the examples presented in Section VI-B. D. Consideration of Uncertainty The proposed method uses load and generation forecasts to dispatch DERs. Therefore, the impact of forecast errors on the accuracy of the proposed method should be evaluated. For this purpose, the following analysis is conducted: loads of Fig. 5 are connected to the grid. For each load a forecast error is assumed that fits a normal distribution with mean value equal to 0 and standard deviation equal to 2.5%. Using these settings, a set of 10.000 simulations is conducted, to statistically quantify the impact of forecasts errors on the accuracy of the proposed method. For each simulation, the power injection at the POI is computed for a wide range of frequencies. The results are depicted in Fig. 8. Specifically, in this figure, the TSO-defined P(f)droop curve is presented (dark blue curve) together with the mean value of the power injections at the POI (black dashed line) and the 95% confidence interval (black ink). As shown, the mean value of power injected at the POI practically coincides with the one required by the TSO. Only trivial differences are observed. The confidence interval remains constant for the entire frequency range. Fig. 8. Impact of forecast errors on the aggregated P(f)droop curve. Fig. 9. Daily load and generation profiles. E. Time Series Simulation In this section, time series simulations are conducted for a typical day, assuming a tequal to 15 minutes. Thus, total a number of 96 time instants is considered. To perform the time series simulations, the load and wind generation daily profiles of Fig. 9 are used. These profiles are scaled to the rated power of loads and DERs of Fig. 5. First, the PFR capability of the examined ADN is evaluated every 15 minutes (t=15 minutes), using Algorithm-1 and the maximum and minimum power that can be provided by the ADN at the POI is computed. Maximum and minimum power limits are presented in Fig. 10. Note that for the examined frequency range, i.e., from 46 Hz to 52 Hz, maximum and minimum power limits are practically constant and not affected by frequency variations. Thus, only two curves are represented in Fig. 10. As shown, power limits vary considerably during the day due to the varying generation and consumption. Therefore, the overall PFR capability of the ADN is also timevarying. The minimum power provided at the POI is always constrained by Pmin,k. The maximum power is generally limited by the available power of the installed DERs, i.e., PMPP,k. However, for certain time periods, e.g., around t=12 h, the maximum power is defined by the thermal limits of the lines. Indeed, Pmax 0is restricted by the proposed Algorithm to avoid congestion issues. To demonstrate Algorithm-2, it is assumed that the TSO requests a P(f)droop curve defined by the following set points: (fmin,0.9×Pmax 0),(fn,0.7×Pmax 0)and (fmax,0.5×Pmax 0). All DERs participate in PFR and wline loss =1. Fig. 11 represents KONTIS et al.: PROVISION OF PFR AS AS FROM ADNS TO TS 4979 Fig. 10. Maximum and minimum power at the POI. Fig. 11. P(f)droop curve required by the TSO and the one provided by DER GWT3. the evolution of the P(f)droop curve settings. In particular, the settings of the TSO-defined droop curve and the settings of GWT3 are presented. As shown, the settings of the TSO- defined droop change continuously during the day and DERs react to these changes to ensure a specific power at the POI. To highlight the time-varying droop settings, droop curves for time instants t1,t2and t3are depicted in Fig. 12. The computational burden of the proposed framework is quantified by calculating the required execution time. The average execution time, required to evaluate the PFR capability of the ADN and to optimally dispatch DERs is 1.64 s and 4.43 s, respectively. This low execution time reveals that the proposed framework can readily support real-time applications in MV networks. Both Algorithms contain steps that can be executed in parallel. For instance, in Algorithm-2, three independent OPF problems are solved. A parallel implementation of these steps will further reduce the computational burden. VII. SIMULATION RESULTS ON A DOUBLE FEEDER DISTRIBUTION SYSTEM In this section, the proposed methodologies are evaluated using a larger distribution system. For this purpose, the 20 kV, 50 Hz MV distribution system of [41] is used. The single line diagram of the examined system is depicted in Fig. 13. As shown, the grid consists of two main feeders. Both feeders, are connected to bus 1. An ITF is used to connect the MV grid to the TS. For the analysis, both PVs and WTs are considered. The rated power of the installed DERs is presented in Table III. The rated active power of system loads is presented in Table IV. All loads are assumed of constant power and they are operated with a power factor equal to 0.9. Further details concerning grid parameters can be found in [41]. The minimum permissible power of WTs and PV units is assumed equal to 10% of Fig. 12. P(f)droop curve required by the TSO and the one provided by DER GWT3 at time instants t1,t2and t3. TABLE III RATED ACTIVE POWER OF PV UNITS AND WTS TABLE IV RATED ACTIVE POWER OF SYSTEM LOADS the Pr,kvalue and zero, respectively. fmax and fmin are assumed equal to 52 Hz and 46 Hz, respectively, while fis set to 0.4 Hz. Vmin and Vmax are considered equal to 0.95 p.u. and 1.05 p.u., respectively. The maximum current of the ITF is 578 A. Thermal limit for the rest of grid lines is 289 A. A. Time Series Analysis To conduct time series analysis, the load and generation profiles of Fig. 9 are scaled to the rated powers of the installed loads and DERs. Initially, Algorithm-1 is executed to determine the PFR capability of the examined ADN. Simulations are performed assuming a tequal to 15 minutes. The maximum and minimum power that the ADN can provide at the POI are illustrated in Fig. 14 together with constraint 8. The