Quantifying the effects of medium voltage–low voltage distribution network constraints and distributed energy resource reactive power capabilities on aggregators

Distributed energy resources (DER), such as, photovoltaic systems and batteries, are becoming common among households. Although the main objective is reducing electricity imports (bills), they could also provide system-level services via an aggregator. However, the more DER provide services, the more important is ensuring that the corresponding operation does not result in network issues. To help DER aggregators understand the implications of network constraints, an AC optimal power ﬂow-based methodology is proposed to quantify the effects that three-phase low voltage (LV) and medium voltage (MV) network constraints can have on the volume of services that can be provided for a given horizon, and the potential beneﬁts from using DER reactive power capabilities. Using a convex multi-period formulation that avoids binary variables for batteries and incorporates voltage-dependent load models, the methodology maximizes DER exports (services) for service-related periods and household self-consumption for other periods (reducing bills). Different service periods are assessed to explore the extent of services throughout the day. Results using a realistic UK MV-LV network with 2400 + households, show that aggregator services can be highly overestimated when neglecting MV-LV network constraints, are inﬂuenced by voltage-demand load characteristics, and that exploiting DER reactive power capabilities can signiﬁcantly unlock further services.


INTRODUCTION
Distributed energy resources (DER), such as photovoltaic (PV) systems, are becoming common in many countries, particularly, in low voltage (LV) networks. This allows households supplying part of their own demand and export the surplus. However, as household exports are typically valued to a lower price than their imports, some are also starting to adopt battery energy storage (BES) systems. This enables households to maximize their selfconsumption by storing the excess PV power and using it when needed. Moreover, given their dispatch capabilities, these households could also participate in providing services to the system. Nonetheless, as these individual DER have relatively small capacity, to effectively participate in the provision of services, they are likely to do so through aggregators [1,2].
The role of DER aggregation to support the system operation has received increasing attention in the last decade, with FIGURE 1 Overview of the proposed methodology load models [25,28], linearizing the corresponding formulation where needed.
The true volume of services (constrained by MV-LV voltage and thermal limits) from a DER aggregator is quantified by maximizing controllable DER exports to provide services for pre-defined service-related periods (hereafter called "service periods") and household self-consumption for other periods (reducing household bills). Different service periods are assessed to explore the volume of services that can be provided throughout the day.
The proposed methodology is demonstrated with a realistic, unbalanced UK MV-LV network with 2400+ single-phase residential customers, where half are considered to have PV and BES systems, and adopting individual (customer level) 30 min resolution time-varying load and PV generation profiles.
This work has been developed as part of the UK EPSRC "MY-STORE" project [29].

METHODOLOGY OVERVIEW
This section introduces the proposed AC OPF-based methodology that allows aggregators to quantify the effects that LV and MV network constraints can have on the volume of services that can be provided (for a given horizon) as well as the extent that DER reactive power capabilities can increase the provision of those services. This quantification is meant to be done offline, so aggregators can decide the most suitable actions ahead. An overview of the methodology and different steps involved are summarised in Figure 1 and described below. For a given horizon, the service periods and selfconsumption periods need to be pre-defined. Their duration need to be aligned with the type of service of interest (e.g. 30 min for some frequency response services in the UK [30]). The service periods to be explored are defined in mathematical terms with the subset T serv ⊆ T , where T is the set of all time periods in the horizon of interest (e.g. T could contain the 48 half-hour time periods in a day and T serv could be a given half-hour period). This means that a service will be provided for all periods t ∈ T serv , and at other times (t ∈ T − T serv ) managed households will maximize self-consumption.
The impact of network constraints and the benefit of reactive power capabilities in the volume of services are quantified by investigating the following three cases: without network constraints (no NC), with network constraints (NC), and with network constraints and DER reactive power capabilities (NC-Q). The no NC case is used to quantify the volume of services that aggregators might be able to achieve from their portfolio when the effects on the network are not considered. NC and NC-Q cases, on the other hand, will show the extent to which those services are affected by network constraints and the role of reactive power to mitigate these limitations.
For each case, different sets of parameters need to be provided to the optimization problem. The capabilities of PV and BES systems are needed to model their operation; this includes kVA rating, power factor range, kWh rating, efficiencies, etc. Such data must be provided to the aggregator by the households participating in providing services. Household-level demand profiles can be obtained from historical smart meter data or using tools such as that in ref. [31]. Solar irradiance profiles are also needed to estimate households' PV generation for the day of interest. This can also come from historical data. Finally, to assess the effect of network constraints it is important to have up-to-date network information such as topology of the network, impedances and tap positions. Therefore, an important assumption in this methodology is that the network-related data is provided by the distribution company to the aggregator. While this information exchange might be challenging in different countries/regions, it is a necessary step towards achieving larger volumes of services without affecting distribution network integrity, similar to what needs to be done between transmission and distribution [32]. The next step is to define the set of constraints to be included in the optimization problem (according to the case) before it is run. Three sets of constraints are defined: household related (including PV and BES balancing equations); network related (e.g. voltage drop and current balancing equations); and DER reactive power related. The full problem formulation (including these sets) is detailed in Section 3.
Finally, after the multi-period optimization problem is solved for each of the three cases, the corresponding volume of services are quantified and compared.

DISTRIBUTED ENERGY RESOURCE AGGREGATOR PROBLEM FORMULATION
This section details the formulation of the DER aggregator problem, which corresponds to a convex multi-period quadratically constrained program. The proposed objective function and the corresponding household constraints (power and energy balances, including DER) are presented followed by the three-phase AC OPF formulation (network model). The proposed formulation is a multi-period extension of the non-linear single-period AC OPF formulation proposed in ref. [28], where non-convex constraints are linearized as in ref. [33], where ZIP load models are incorporated with its corresponding convex approximations, and, crucially, where binary variables are not needed to avoid the simultaneous charging and discharging of BES. Thus, improving scalability and solving speed for large unbalanced distribution networks.
The proposed formulation focuses on export services, specifically the injection of active power from managed households to the network; which is aligned with what is being trialled in recent projects due to the nature of existing markets [34][35][36]. Nonetheless, the formulation can be modified for import services (absorption of active power), if desired by the aggregator due to a market need. These modifications are described in Section 3.2.
BES powers can be considered continuous variables. Discrete network components managed by the distribution company, such as on-load tap changers (OLTCs), are assumed to be shared with the aggregator (as previously discussed) and, therefore, are taken as input parameters. Consequently, the optimization problem is only composed of continuous variables. For the benefit of the reader, inputs parameters are shown in bold.

Objective function
The proposed objective function in Equation (1) has three terms and aims to find the largest available aggregated power for pre-defined service periods (t ∈ T serv ⊆ T , where T is the considered horizon, e.g. a day), while minimizing household imports (bills) for the rest of the periods (t ∈ (T − T serv )).
. (1) The first term in the objective function is implemented by minimizing the net household power consumption (imports minus exports), P s,t , over all households managed by the aggregator (set S ⊆ H , indexed by s, where H is the set of all households, indexed by h) and for times t ∈ T serv . This intrinsically minimizes household power imports, P + s,t , and makes power exports, P − s,t , as large as possible (see Equations (5) and (7)). The second term is implemented by minimizing household power imports, P + s,t , and exports, P − s,t , over all households for times t ∈ (T − T serv ). This results in a minimization of network dependency for managed households, that is, maximization of households' self-consumption. Weighting factors ( serv ≥0 and self ≥0) can also be considered to give priority to the volume of services during T serv ( serv > self ) or to minimize household imports for the rest of the day ( serv < self ). This prioritization is relevant in cases where household resources are not enough to fully avoid imports. Moreover, these weighting factors can have different values in time to give priority to certain service periods ( serv t ) but also to model the actual selfconsumption operation of the BES systems ( self t ) when not providing services; that is, charge with surplus generation, discharge to supply demand [37]. The latter is implemented with a gradually decreasing weighting factor ( self t > self t+1 ) to ensure that as soon as demand can be supplied or generation stored, the BES system operates accordingly.
The third term is introduced to avoid solutions where BES systems (set B, indexed by b) can have charging and discharging powers simultaneously (P b,+ s,t and P b,− s,t , respectively). With > 0, the non-negativity of these two variables (see Equation (7)) makes their sum to be minimal when at least one of them is zero. For instance, a net charge of X kW (see Equation (6)) for the managed household s at time t , can be achieved with and for any ≥0 (where simultaneous charging and discharging occurs for > 0). Since any > 0 needs to be avoided, the third term in the objective function, (X + 2 ), is built such that it forces = 0 when minimized (note that the same applies to a net discharge of X kW). This effectively avoids the simultaneous charging and discharging of BES systems without the need for binary variables, which are commonly used in the literature (e.g. [7,11,38]). This, in turn, results in a formulation made only of continuous variables, that is, more scalable and much faster to solve. The penalization factor should be small enough to give priority to the first two terms, but large enough to avoid making the third term irrelevant (where simultaneous charging and discharging could occur). This point is further demonstrated in the Section 5 for different values of . Note that the final convex problem (shown latter), can also be written as a maximization problem by making the objective function to be concave (multiplying Equation (1) by −1).

Household and distributed energy resource constraints
The active power balance of each managed household at time t , P h s,t , is given by Equation (4), where P d s,t and P g s,t are the household non-controllable (by the aggregator) demand and PV generation, respectively, and P b s,t is the total power of each controlled BES. The fact that P d h,t could vary with the voltage magnitude has been represented in Section 3.3.2, through the use of ZIP models.
As previously mentioned and to facilitate the formulation of the objective function, P s,t , which can take positive or negative values, is split into imports (P + s,t ) and exports (P − s,t ) in Equation (5) with non-negative variables Equation (7). Similarly, P b s,t is split in Equation (6) into charging (P b,+ s,t ) and discharging (P b,− s,t ) non-negative variables (see Equation (7)). Another important constraint that makes sure every managed household gets to provide services (injecting power) is Equation (8). This ensures that all these households will either discharge their BES or stay idle at service periods. Note that for import services, only two changes are required in the formulation: the first term in the objective function and Equation (8) have to be multiplied by −1.
The inter-temporal constraint for the calculation of the energy stored in a BES system, E b s,t , which depends on its previous charging/discharging powers as in Equation (9), considers charging and discharging efficiencies, b,+ s and b,− s , respectively. Here, E b s,0 is the energy stored at the BES system at the beginning of the day and Δt is the time-step used in the multi-period optimization. Also, the energy stored in the BES system is limited by its energy rating, E b s , as in Equation (10), and its depth of discharge, DOD b s , as in Equation (11).
Finally, the active and reactive power output (capabilities) of BES systems, P b s,t and Q b s,t , respectively, are limited by the apparent power rating of their inverters, S b s , and it is given by Equation (12). It should be noted that constraint Equation (12) assumes that the BES system has a dedicated inverter. For hybrid PV-BES systems where a single inverter is used, Equation (12) should be substituted with Equation (13), where Q inv s,t and S inv s are the reactive power output and apparent power rating of the hybrid inverter.
( P g s,t + P b

Three-phase AC optimal power flow formulation
This section presents the linear and non-linear constraints of the original AC OPF formulation in ref. [28], as well as, the convex approximations, as in ref. [33] for shunt currents and minimum voltage limit, where its accuracy was demonstrated. Additional convex approximations are proposed for incorporating ZIP load models to the problem. The resulting formulation allows to model the effect of network constraints (voltage and thermal limits) for large MV-LV networks in the aggregator problem. For brevity, the time-dependent subscripts, which apply to currents, powers and voltages, are omitted.

3.3.1
Linear constraints The following constraints are linear (convex) in the original formulation [28] and, therefore, do not need approximations. Constraints in Equations (14) and (15)  are the resistance and reactance of the line between phases , ∈ Φ. Shunt admittances have been neglected due to their typically negligible values [28].
(15) The current balance at each node n ∈ N and phase ∈ Φ, given by Kirchhoff 's Current Law and shown in Equation (16), is also a linear constraint. This considers all shunt and series elements connected and is directly split into real and imaginary parts. Here, e indexes lines and transformers in sets L and Y , respectively. Indexes o, g, b, and d represent the source generator (transmission interface), distributed generators, BES systems and loads, respectively. Here, BES currents take positive and negative values when charging and discharging, respectively.
For each delta-wye 11 (Dy11) transformer (y ∈ Y ), typical in distribution networks [39], the equivalent phase to ground complex voltage relation is given by Equation (17), and the phase currents are related by Equation (18), as deduced in ref. [39]. Here, 1 y and 2 y are the delta and wye nodes, respectively, + 1 and + 2 are the two consecutive phases to (e.g. for = 2, these are 3 and 1, respectively), Z y is the complex impedance of the transformer (Z y = R y + jX y ), r y is the transformer nominal ratio (equal to √ 3 if the voltage at each side coincide with the network nominals), and y is the three-phase transformer tap position in pu.
In the original problem of ref. [28] (managed by the distribution company), y were variables to be optimized (for OLTCs). However, these cannot be controlled by the aggregator who takes them as input parameters. This makes Equations (17) and (18) to be linear (convex) for this problem. Finally, Equations (17) and (18) are split in real and imaginary parts.

Non-linear constraints and their approximations
These correspond to the representation of voltage dependent load powers, shunt element currents components and operational limits.
• Voltage dependent load powers: The voltage dependent active power, P d h , and reactive power, Q d h , demanded by a single-phase load d ∈ D at household h ∈ H (managed by the aggregator or not), connected at node n and phase , are represented in pu using ZIP models in Equations (20) and (21). Here, P d0 h and Q d0 h are the active and reactive power at nominal voltage and Finally,Ẑ are the constant impedance, constant current and constant power parameters for the active or reactive power components [25,28].
In the objective function, P d s (for managed households only) is non-convex due to V ) as shown in Equation (22), and replaced back in Equation (20).
• Shunt element currents: The instant real and imaginary parts of the current demanded by a single-phase household load d ∈ D at household h ∈ H (managed by the aggregator or not), connected at node n and phase , are represented in pu by Equations (23) and (24).
where P d and Q d are given by in Equations (20) and (21). The resulting current equations, are linearized with a first order Taylor expansion around (V re n, * , V im n, * ) as in Equation (25).
Single-phase PV systems g ∈ G are assumed to operate at unity power factor, they can be approximated as constant power devices [25], and are not managed (i.e. not curtailed). Thus, time-varying PV powers are considered to be "known values" (e.g. based on forecast). Their currents can be approximated with a similar set of linearized equations. Note that the power of non-controllable elements (e.g. constant power loads and PV) can be aggregated to obtain their corresponding currents, reducing the number of variables of the problem.
BES systems (b ∈ B) are assumed to be directly managed by the aggregator and for which the benefits of regulating reactive power are also explored. The active and reactive powers for a BES connected at node n and phase are approximated using Equations (26) and (27), respectively (positive when charging and with Q b s = 0 when not managing reactive power).
Note that three-phase connected customers, with or without embedded generation, at MV or LV networks, can be represented with the same set of equations per phase.
The accuracy of these linear approximations will depend on the chosen operating point, which can be estimated using recent operational data. To further improve accuracy (if needed), a second multi-period OPF could be run, where the voltage results of the first serve as proxy operating points for the second [33,40]. Another simple approach can be using expert knowledge to choose voltage operating points. For instance, household voltage magnitudes closer to the upper limit at service periods (maximum services), and around the transformer secondary side voltage (determined by the tap position) at self-consumption times (due to low power flows).
• Network limits: The magnitude of phase currents in each line l (set L) and transformer y (set Y ) are constrained by their capacity I e , with e ∈ L ∪ Y . These constraints can be expressed in a quadratic and convex way by relating their square magnitudes as shown in Equation (28).
where 1 e , 2 e ∈ N are the starting and ending nodes of e. The constraint for the maximum voltage limit at node n and phase , V n , is represented similarly with Equation (29).
However, the analogous constraint for the minimum voltage limit, V n , is non-convex. Here, the voltage magnitude squared ) is linearized with a first order Taylor expansion around (V re n, * , V im n, * ), resulting in the linear constraint shown in Equation (30).
As in all previous linearizations, the accuracy of Equation (30) depends on the quality of the estimated voltage. However, the maximum voltage limit constraint, likely to be the most enforced when providing services, is not linearized.

CASE STUDY
The proposed AC OPF-based methodology is demonstrated with a realistic UK MV-LV network with 2430 single-phase residential customers and considering that 50% of those customers have PV and BES systems. This large PV penetration does not create thermal or voltage issues, that is, reverse power flows can be handled without the need for new assets or operational strategies. However, network issues can arise if an aggregator discharges BES systems at times of PV generation. The optimization program is implemented in AIMMS with Gurobi 8.1 as solver [41]. Voltage limits in the AC OPF are defined based on the UK standard BS EN50160 [42], which defines that the 10 min average rms customer voltages must remain between 0.94 and 1.1 pu (base 230 V line to neutral). These voltage limits are enforced at all periods. Finally, for the

Realistic UK medium voltage-low voltage network
A residential MV-LV network from the North West of England, shown in Figure 2, is used. The MV and all LV feeders have been validated in close collaboration with the corresponding distribution company, Electricity North West Limited, using monitoring data and according to the criteria described in ref. [43]. This network consists on a 6.6 kV MV feeder that supplies 2430 single-phase residential customers through 9 Delta-wye 11 (Dy11) secondary substations and 36 underground threephase four-wire 0.4 kV LV feeders (transformer capacities and number of feeders per LV network are also shown in Figure 2). This results in a network model of 4579 three-phase nodes, suitable for assessing the implementability of the proposed approach. The main characteristics of the LV networks are summarized in Table 1. For illustration purposes, the topologies of two LV feeders with the smallest and largest number of customers are shown in Figure 3(a),(b), respectively; where the secondary transformer and customers per phase can be seen.

Load, photovoltaic and battery energy storage modelling: Time-varying profiles
A pool of one thousand 30 min resolution domestic household profiles are created using a tool developed in ref. [31], hence, the time-step used in the multi-period optimization is t = 0.5 h. This tool uses UK statistics to determine the status of domestic appliances in a day considering occupancy of the house, type of day and seasonality. To approximate the occupancy of houses when creating these profiles, UK statistics are also used [44]. Then, active profiles are obtained per household by aggregating individual appliances. Reactive power profiles are calculated considering a power factor of 0.98 (inductive). These values are used as the powers at nominal voltage for the ZIP model (P d 0  h,t ,Q d 0  h,t ). These profiles are randomly assigned among households. For the modelled MV feeder (with 2430 households) this results in a peak demand of approximately 2 MVA. For illustration purposes, Figure 4 shows: (a) a single household active/reactive power profile and (b) the diversified profiles considering all one thousand pair of active/reactive profiles.
PV generation is simulated with unity power factor and a linear relationship between PV system ratings and solar irradiance [25]. In this study, 50% of the households are randomly selected to have PV and BES systems. PV systems are allocated with ratings aligned with UK statistics [45] (from 1 to 4 kWp and with an average of 3 kWp). BES systems are considered to be 3.6 kW, 4.2 kWh, with 90% depth of discharge (i.e. SOC lower limit of 10%) and 90% round trip efficiency [46]. A 30-min resolution clear-sky solar irradiance recording from a summer day in Manchester, UK [47] shown in Figure 4(c), is used. For simplicity and given the relatively small geographical distance (less than 5 km), it is assumed that all customers get the same solar irradiance.

Objective function weighting factors
The objective function is implemented with the weighting factors shown in Equations (31)-(32) and = 0.1. These are meant to prioritize customer self-consumption, that is, the aggregator manages only spare capacity.

Quantifying the volume of services
When considering network constraints, the available headroom to provide services will vary throughout the day due to the changes in demand and PV generation. To understand this effect, three service periods of interest are studied, each with a duration of 30 min and first considering constant power load models (the impact of different ZIP load models is investigated in Section 4.6). To explore an early morning service (low demand and low PV generation), the service period (T serv in Section 2) is set to start at t = 13 (6:30 AM), that is, the service is provided between 6:30 and 7:00 h. Then, T serv is set to start at t = 26 (1:00 PM) to explore a service at the time of peak PV generation. Finally, T serv is set to start at t = 42 (9:00 PM) to explore a service at summer peak demand time. For each of these service periods, the volume of services that can be provided are quantified for the three cases defined in Section 2: a) without considering network constraints, "no NC," b) considering network constraints, "NC," and c) considering network constraints and reactive power capabilities from BES inverters, "NC-Q." For comparison purposes, the baseline self-consumption only (SC) operation of BES systems, that is, no services, is also included. To adopt a realistic SOC at the beginning of the day of interest, a prior similar day is also considered in which the BES systems start from the SOC lower limit (10%). The half-hourly average SOC for the second day (all BES systems) is illustrated in Figure 5 which shows initial and final SOC of ≈40%, and that full capacity is reached around 2:00 PM. For consistency, the three cases to be investigated consider all BES systems to have an initial SOC of 40%.
Half-hourly bus voltage magnitudes and line/transformer utilization levels (percentage of phase current with respect to rated capacity) for the SC case are shown in Figure 6(a),(b), respectively. These values help understanding how much   voltages or currents could be increased at different times of the day without violating network constraints, that is, headroom to provide services. The network has the lowest voltage and thermal headroom from 10:00 AM to 4:00 PM due to high PV generation, as expected. Furthermore, the most constraining factor seems to be voltage levels in LV networks, followed by the capacity of secondary transformers (TX).

Services during early morning
The aggregated power (positive for imports and negative values for exports) of all managed households providing services at 6:30 AM is shown in Figure 7 for the cases "no NC," "NC" and "NC-Q." It can be seen that for all of them, the aggregated operation is identical (3 MW of aggregated exports). This means that, at this time, network limits are not reached, hence, the volume of exports that the aggregator can provide are only limited by the stored capacity of BES systems. The latter is demonstrated in Figure 8(a) where the average SOC among all BES systems is shown. They reach the SOC lower limit (10%) during the service period, hence, there is not enough energy stored to provide more active power. Furthermore, it should be noted that the aggregated power (exports) do not necessarily indicate the volume of services that the aggregator could provide. This is because the volume of Aggregated managed households' power per case for T serv = 1:00 PM services is the difference in aggregated power between the case providing services and what would otherwise occur without services, that is, the SC case. This aggregated difference is plotted in Figure 8(b) from 4:00 AM to 12:00 PM, with the service period being highlighted. During the service period, the difference (the volume of services) matches the exports seen in Figure 7, that is, 3 MW. This is because the aggregated power with the SC case was 0 MW (i.e. all managed households are fully self-sufficient). Interestingly, from 8:00 AM to 12:00 PM, a "rebound" effect is seen. Because the BES systems were discharged by 7:00 AM (end of service period), they charge with surplus PV generation as soon as possible to ensure demand is covered later in the day. This, in the aggregate, when compared to the SC case, results in less exports from 8:00 AM to 12:00 PM (Figure 7). While this has no effect on the level of services that can be provided at 6:30 AM, system operators need to understand this effect as it can affect estimates related to PV generation.

Services during peak generation
The aggregated power of all managed households providing services at 1:00 PM for all the cases is shown in Figure 9. At this time, the level of services that can be provided is limited by network constraints as anticipated from Figure 6. The "no NC" case results in a large overestimation of the active power that the aggregator could export (7.1 MW) compared to the value obtained when network constraints are catered for (3.6 MW). Although the latter value is much smaller, it is also demonstrated that aggregated exports can be significantly improved by enabling reactive power capabilities from BES systems (to 4.8 MW, a 33% increase). Note that the benefit from enabling BES reactive power capabilities become even clearer when these aggregated exports are translated into volume of services. At this time, the aggregated power for the SC case is 2.7 MW of exports. Then, the volume of services for the "NC" and "NC-Q" cases are 0.9 and 2.1 MW, respectively. Therefore, when considering network constraints, the use of reactive power capabilities from BES systems can result in more than the double of the volume of services compared to the "NC" case, as it can be seen for this service period in Figure 11(b). The significant difference between the NC and NC-Q cases is attributed to how reactive power capabilities enable a better utilization of the distribution network infrastructure. This is demonstrated in Figure 10 where the different variables that can constrain the network (voltages and line/transformer utilization levels) are shown for the service period. In both the NC and NC-Q cases there are node voltages, as well as some lines and transformers, at their maximum limit. However, in the NC-Q case, this occurs to many more nodes, lines and particularly transformers mainly due to the larger amounts of active power exports being achieved. In the NC-Q case, reactive power is absorbed or injected by the BES systems (as shown in Figure 11(a)) where needed. In areas where voltage limits are binding, reactive power will be absorbed to counteract the voltage rise caused by active power exports. In areas where thermal limits are binding, reactive power will be injected by other BES systems to counteract reactive power flows (absorption from demand or BES systems), making the most of thermal capacity for active power flows. Both of these effects are seen at 1:00 PM (see Figure 11(a)) as voltage and thermal limits are both binding; 53% of BES systems end up absorbing reactive power (1.6 Mvar) and 47% injecting reactive power (1.3 Mvar). This demonstrates that by enabling BES reactive power capabilities it is possible to unlock the provision of more active power services.
As in the early morning service, the "rebound" effect is also evident for this service period, as shown in Figure 11(b), where bars overlap. However, it is different in its extent and when it happens. The larger the volume of services, the more dramatic battery energy storage systems' active power for T serv = 9:00 PM this rebound effect is. Also, this effect happens right after the service is provided due to the available surplus PV generation that allows charging the BES systems (SC case in Figure 9).

4.4.3
Services during peak demand Figure 12(a) shows that for the peak demand service period (i.e. 9:00 PM), the aggregated managed household' active power that could be provided as service is the same for all cases (2.7 MW). This is because network limits are not reached, similar to the morning service period. In addition, since in the SC case the households manage to fully supply their demand from BES systems at 9:00 PM (zero net demand), the volume of services also corresponds to 2.7 MW. The volume of services during this period is limited by two aspects: the available headroom of BES inverters and the spare stored energy. Because BES systems need to discharge to cover the demand, when required to provide services, exports are limited by the available headroom (full BES inverter capacity minus demand being supplied). In addition, exports will also be limited to the spare energy (energy to cover the subsequent demand period minus stored energy). Figure 12(b) shows the active power of the BES systems for the SC and no NC/NC/NC-Q cases (same for all). The SC operation represents the demand that needs to be met (median of 0.45 kW, consistent with Figure 4(b)). When providing services, half of the households with enough spare energy use their full BES inverter capacity to export power. However, those houses with small volumes of spare energy end up exporting less power.

Single-period availability of services across the day
This section extends the previous analyses by assessing the maximum available volume of services that could be provided by the aggregator for a single time period but with results presented throughout the day. This is obtained by applying the proposed OPF-based methodology to each of the time steps in the case study (i.e. 48), where T serv takes only one of the time steps per simulation. By obtaining the maximum volume of services at each time step (aggregated managed household power minus their SC power) it is possible to build "envelopes" for each case, as shown in Figure 13. Here, constant power load models are used for consistency with Section 4.4.

FIGURE 13
Single-period availability of services across the day, per case During the early hours of the day (before 5:00 AM), the maximum volume of services that can be provided by the aggregator is relatively constant (1.3 MW of exports) and the same for all cases. This is because at this time they are only limited by the energy stored in the BES systems (average SOC shown in Figure 5) as charging from the grid is penalized in the objective function, hence, charging occurs with surplus PV generation only. Between 5:00 AM and 6:00 AM, when PV systems start generating, the maximum volume of services increases as BES systems start charging, leading to more energy stored. Between 6:00 AM and 7:00 PM, a significant difference on the maximum volume of services can be observed between the three cases, with a significant overestimation when not considering network constraints. This highlights that network issues are likely to occur when providing services without considering network limits. This is due to the PV generation causing higher voltages and line/transformer utilization, leaving a smaller headroom for the BES systems to inject active power. This effect reaches its peak at around 1:00 PM, as expected due to peak PV generation. During high PV generation, reactive power capabilities from BES systems can significantly unlock services, with a maximum increase of 33% at 1:00 PM, as shown in Section 4.4.2.
After 7:00 PM (peak demand and little PV generation), all three cases result in the same available volume of services (around 2.7 MW) as network constraints are not binding. During the evening, it is the need for BES systems to cover demand that constrains their service provision.
Finally, the largest service availability occurs from 7:00 AM to 9:00 AM, as BES systems make the most of surplus PV generation (see Figure 5), resulting for many households in zero net demand. Given that upcoming periods also have surplus PV generation (that can be used to charge again), BES systems are in a position to go from charging to discharging. This, in effect, provides the largest possible volume of services as it combines the injections from PV systems and BES systems.

Load model impact
This section explores how different voltage dependent characteristics of loads can affect this potential volume of services. The three main load models are considered: Constant power (as in previous sections), constant current and constant impedance (their weighted combination is also possible). For simplicity, active and reactive load powers are assigned the same load model and for all households (managed by the aggregator or  ). Therefore, as all loads can be affected by voltages, the quantification of the potential volume of services is assessed by adding up all households' net power for the NC (or NC-Q) case, and subtracting the corresponding SC case value. The case without considering network constraints has been omitted as voltages largely outside the statutory limits are not acceptable and loads would be unrealistically affected.
Here, the peak generation service time (same as in Section 4.4.2) has been selected as this was shown to be the most affected by network constraints (the volume of services differed the most among cases). Table 2 shows the sum of all households' power (MW) and the corresponding volume of services that could be provided when considering each of the three load models, for the NC and NC-Q cases. The average voltage among LV customers (in pu) is also included for reference.
For each case, the resulting voltages are, in average, higher than the nominal and with similar values regardless the load model. This is because power flows (reverse) are dominated by DER powers at this time (constant power devices). Moreover, voltages are higher when providing services (NC) and even higher when using reactive power capabilities (NC-Q), due to larger reverse flows, as expected.
For the SC case, voltages are already high due to PV-induced reverse power flows. This makes loads with larger voltage dependency consuming more (constant Z and constant I cases) and, therefore, reducing the network reverse power flow (exports due to PV only). On the other hand, when providing services, the larger the DER exports, the higher the voltages. Again, this makes loads with larger voltage dependency consuming more, which brings voltages down allowing some extra room for more DER injections (optimized) across the network. This translates to larger reverse all households' power and larger volume of services.
This effect is more evident for the NC case, where the increase between the volumes of services when considering constant impedance loads compared to the case of constant power is 44%. For the NC-Q case, although the volume of ser-vices is larger for any load model (compared to the NC case), the same difference represents only 7%. The latter is explained for voltages being already optimized (with reactive power compensation) for all load model cases.
Note that, for the NC and NC-Q cases, the volume of services difference among load model cases is larger than the corresponding all households' power difference. This is because for the SC case, the all households' power decreases (less reverse powers) from the constant power to the constant impedance case, as previously explained, but it increases in the same direction for the NC and NC-Q cases.

DISCUSSION
In this work we proposed a methodology to estimate the volume of services that can be provided for a single aggregator. This offline quantification could also be applied when multiple aggregators are in place, assuming that each will try to maximise their volume of services and, therefore, they all can be treated as one [20]. However, for this assumption to hold there are some challenges, such as, each aggregator having the same projected information of the DER managed by others (load and generation forecasting). The interactions of different aggregators during actual operation (online) is beyond the scope of this work.
Although the methodology is for offline purposes, it can provide the basis for an extended version for operational (online) purposes. Aspects of demand/generation forecasting, communications and DER dispatch fairness should be considered among with modelling the interactions between different aggregators. To cater for load and PV generation uncertainties, the methodology can be applied, for instance, under many different random (Monte Carlo approach) or extreme load/PV generation scenarios, informing the aggregator on the volume of services it could offer under different conditions. The aggregator can adopt a more or less conservative approach, however, the methodology would be unaffected on its formulation.
All the network impedances are considered to be known a priori by the distribution company and shared with the aggregator as this benefits both parties. Given that, in reality, actual impedances can present variations with respect to informed values, a sensitivity analysis could be carried (e.g. as in ref. [48]). Similarly, the topology of the network is assumed to be known and fixed throughout the analysis. However, the model can be updated to cater for network reconfigurations or expansions, provided the corresponding data is made available to the aggregator.
If the network to be modelled contains capacitor banks, these can be modelled as negative constant impedance loads [28], non-controllable by the aggregator.
Although the proposed objective function can be classified as a weighted sum approach, it does not correspond to three separate and conflicting objectives. In fact, the first and second terms operate over different time sets and are closely related, while the third term aims to avoid the simultaneous charging and discharging of BES systems at all periods. However, if the energy exports (first term), self-consumption (second term) and penalization of simultaneous BES charging and discharging (third term) are "converted" into, for instance, costs, there might be an opportunity to explore potential trade-offs (Pareto approach). Finally, it is demonstrated with some case examples that different penalization factors, , can produce the same performance for aggregators while avoiding simultaneous charging and discharging of BES systems, as far as they are not too small to make the third term irrelevant. Table 3 shows the resulting aggregated managed households' power (AP (MW)) for equal to 0.9, 0.1 and 0.01, and for the SC, NC and NC-Q cases. The performance of the third term avoiding simultaneous charging and discharging of BES systems is indicated with the average of the multiplication in pu of these two variables across the day and among all BES, as in Equation (33).
For simplicity, all loads are assumed to be constant power type. It can be seen that while the performance for the aggregator is nearly the same in all cases, the smallest (0.01) is too small already to guarantee avoiding simultaneous charge and discharge for BES. Therefore, a sensitivity analysis is suggested for choosing .

CONCLUSIONS
This work proposes an AC OPF-based methodology to help DER aggregators quantify and understand the effects that three-phase LV and MV network constraints can have on the volume of services that could be provided from their DER portfolio throughout the day. The methodology also allows quantifying the benefits from using DER reactive power capabilities and the impact of different loads' voltage dependency (load models). A convex multi-period formulation (suitable for large networks) is used to maximize DER exports for service periods and household self-consumption for other periods. In this formulation, solutions with simultaneous charging and discharging of BES are avoided without the need of binary variables. Results in a realistic UK MV-LV network with 2400+ households, half of them having a PV+BES system, demonstrate that aggregator services can be highly overestimated when neglecting MV-LV constraints, particularly at times of high PV generation. This is because reverse power flows due to PV generation leave little room for further exports from BES systems before reaching network limits. However, further volumes of services can be unlocked by enabling reactive power capabilities from BES systems. In areas where voltage limits are binding, reactive power will be absorbed to counter-act the voltage rise caused by active power exports. Where thermal limits are binding, reactive power will be injected to counteract reactive power flows (to demand or BES systems) making the most of thermal capacity for active power flows. At times of little PV or no generation, the volume of services is not constrained by network limits but by the energy and power capabilities of BES systems; the former due to the evening demand that needs to be covered and the latter due to the headroom left from supplying the local demand. The largest headroom to provide services (without reaching network limits) occurs at times when BES systems charge from surplus PV generation (exports from PV and BES systems combined).
Results also show that the loads' voltage dependency can have a large impact on the quantification of the volume of services when not using DER reactive power capabilities. This is because, for voltage dependent loads, an increase in voltages due to reverse power flows, lead to higher demand which allows for some extra room for more DER injections (optimized) across the network. However, this effect is less important when enabling DER reactive power capabilities given that, in this case, voltages are being already optimized (with reactive power compensation) for all load model cases.
Finally, for services provided during PV generation, a rebound effect occurs as BES systems charge again from surplus PV generation, hence, increasing net demand from aggregator-managed customers. System operators will need to understand and anticipate this effect following service requests from this type of aggregators.  Real/imaginary current for element e at connection node e ∈ N and phase ∈ Φ V re n, , V im n, Phase to ground real/imaginary voltage at node n ∈ N and phase ∈ Φ