Exploiting OLTC and BESS Operation Coordinated with Active Network Management in LV Networks

Abstract

The large number of small scale Distributed Energy Resources (DER) such as Electric Vehicles (EVs), rooftop photovoltaic installations and Battery Energy Storage Systems (BESS), installed along distribution networks, poses several challenges related to power quality, efficiency, and reliability. Concurrently, the connection of DER may provide substantial flexibility to the operation of distribution grids and market players such as aggregators. This paper proposes an optimization framework for the energy management and scheduling of operation for Low Voltage (LV) networks assuring both admissible voltage magnitudes and minimized line congestion and voltage unbalances. The proposed tool allows the utilization and coordination of On-Load Tap Changer (OLTC) distribution transformers, BESS, and flexibilities provided by DER. The methodology is framed with a multi-objective three phase unbalanced multi-period AC Optimal Power Flow (MACOPF) solved as a nonlinear optimization problem. The performance of the resulting control scheme is validated on a LV distribution network through multiple case scenarios with high microgeneration and EV integration. The usefulness of the proposed scheme is additionally demonstrated by deriving the most efficient placement and sizing BESS solution based on yearly synthetic load and generation data-set. A techno-economical analysis is also conducted to identify optimal coordination among assets and DER for several objectives.

1. Introduction

Low Voltage (LV) distribution networks used to be a passive segment of the power system, mainly for the supply of consumers; thus, power flows were heading from the bulk transmission points to the distribution grid. Accordingly, from the secondary substation and the downstream connected LV grid, there used to be very limited or an absence of automation for its monitoring and control [1]. In the last decade, there has been a large number of small-scale units, commonly referred to as Distributed Energy Resources (DER), that are getting connected along distribution grids. Several types of DER may be connected such as domestic rooftop Photovoltaics (PV) or generally microgeneration (μ G) in some cases coupled with Battery Storage Systems (BESS), controllable loads (e.g., Electric Heat Pumps or other smart domestic appliances), and Electric Vehicles (EVs). The extensive integration of DER in the grid may cause several technical challenges on the operation of distribution networks such as voltage problems, branch congestions, and phase unbalances. Despite these technical bottlenecks, DER may be utilized in favor of the grid operation, providing ancillary services and supporting the bulk transmission system and distribution networks [2,3].

The Distribution System Operators (DSOs) are currently adopting practices to enhance the observability and controllability of the distribution grids throughout Advanced Distribution Management Systems (A-DMS), [4]. The active involvement of DER and generally prosumers in the operation of the network is generally referred to active network management, which is regarded to take place utilizing their flexibility. Sources of flexibility may come from several types of DERs that are enabled with temporal shifting of active or reactive power to be consumed or injected into the grid. Such strategies of active participation of consumers in the grid’s operation have gained the interest of utilities for the past few decades by engaging, mainly, industrial consumers through demand side management schemes [5]. Several research works have discussed recently about the Smart Transformers (STs) envisioned as a key element for the controllability of distribution networks in a future context of DER massification [6]. For the smart grid development, more novel advanced control schemes have to be implemented towards the active involvement of DER.

2. Related Works and Contributions

The operational control schemes and energy management applications could be generally classified into several categories according to the communication infrastructure and the data requirement (i.e., deemed necessary to be used). Based on the latter, centralized schemes usually look for solutions not only to resolve technical grid constraints but also to optimize the economical operation of the grid [7,8,9]; local control (or decentralized) techniques may be applied merely relying on droop based rules [10,11] and distributed strategies which are in line with the deployment of local energy communities and transactive energy concepts [12,13,14].

In LV distribution networks, voltage regulation and phase balancing are managed by the DSOs, typically by manual adjustments (offline) of the MV/LV—secondary—transformer which may happen once or twice a year, depending on the seasonal changes in the loads [15]. Alternatively, DSOs act by investing on grid reinforcement measures such as line replacement (i.e., when branch congestions) and manual phase redistribution for phase balancing [16,17]. Considering the stochasticity of both load and generation, the aforementioned practices of manual configuration of tap-positions and grid reconstruction may be inadequate in many cases [18]. Manual controls and simple local controls may be insufficient due to the intermittent nature ofμG and the stochastic behavior of EV charging. On the other hand, the grid reinforcement may be considered quite effective but still a costly measure for the DSO.

The possibility of utilizing droop capabilities (for active and reactive power control –P=f(V),Q=f(V) , accordingly) with a smart PV inverter particularly for voltage regulation has extensively been studied in the literature [10,19,20]. The reactive power control is generally a less efficient solution in the LV grid for voltage control due to the high branch ratio R/X (i.e., rather resistive nature of LV distribution lines) compared to Medium Voltage (MV) distribution networks or transmission. Self-consumption is commonly imposed by regulation and legislation lately, to address voltage rise effects during the peak period of PV generation. In several European countries (e.g., Belgium, Denmark, the Netherlands and Greece), residential PV self-consumption measures based on net metering schemes aim at matching the endogenously generated power with local demand [21]. In Germany, there is a cap for active power feed-in at 70% of the installed capacity for all the prosumers with a capacity of less than 30 kWp [22]. Nevertheless, Active Power Curtailment (APC) might not be an economically attractive solution for both DSOs and the prosumers. Therefore, more sophisticated control schemes are proposed exploiting the coordination ofμ G with DSO assets to improve the network’s power quality [23,24].

Several applications have focused their interest on introducing the control of other DER, such as BESS or EVs. Most research applications refer to the coupling of BESS systems withμ G to firm-up the dispatched power produced from PVs by reducing the mismatches between generation-demand [21,25]. The increased cost of investment has been the main limitation for the extensive deployment of BESS, a fact that is likely expected to change in the current decade according to [26]. Concurrently, several works have lately proposed the utilization of BESS by DSOs—i.e., owned and controlled by the DSO—[27,28] to deliver operational flexibility operation as well as to increase hosting capacity of DER into the grids. Nonetheless, there is very limited BESS utilization by DSO currently due to the in force Directive 2009/72/EC [29], where the unbundling requirements for DSOs do not allow BESS directly owned and controlled by them. As a result, the ever growing number of domestic BESS may undermine the current business model of the electric utilities [30]. The following trend aims to maximize the revenue brought by “smart” consumers that utilize home energy management systems to optimize the local generation and consumption.

The large-scale penetration of EVs that is expected in the current decade will notably increase electricity consumption, during charging periods. Therefore, power flows—including Vehicle-to-Grid (V2G)—grid losses, and voltage profile patterns and generally power quality along the grid will change significantly [31]. These effects may arise the need to reinforce the grid in some locations. Based on the EV charging strategy to be adopted, grid reinforcement may be deterred. Several schemes have been proposed to derive smart charging schedules to ensure safe grid operation [32,33,34,35], while some of those dealing with phase unbalances may be provoked by EVs in the grid. None of the aforementioned works, however, propose any possible coordination of DER amongst them or with DSO assets to optimize cost objectives or the technical operation.

Several industrial prototypes for secondary transformers (Efacec, Reinhausen, Siemens) are equipped with the capability of On-Load Tap Changer (OLTC) [36] for MV/LV transformers. There is relatively limited work dealing with the coordinated operation of OLTC with DER addressed in [30,37,38,39]. In spite of the fact that these works provide the optimal coordination of OLTC with DER, there is no insight for the temporal flexibility that may be delivered by DER between subsequent time slots. Authors in [7] propose a framework for the optimal coordination among several DER and the OLTC, dealing also with the phase balancing constraint. In this work, authors propose efficient linearizations to resort tractable multi-period OPF extending the problem statement in [40]. On the contrary, a three-phase multi-period OPF based on the exact (i.e., nonlinear) AC power flows is proposed in this work, incorporating multiple DER within the operation of the distribution grid.

This paper advances previous works of the authors [41,42]. The main contributions of this paper may be outlined as follows:

1. Advances an analytical DMS framework for the energy management and scheduling of operation of unbalanced distribution networks with increased integration of DERs. The tool is capable of deriving control actions and schedules for flexible DER and the OLTC subjected to multiple operational constraints such as congestion management, phase balancing, and voltage regulation. Furthermore, optional objective terms might be opted among the minimization of operational costs or minimization of flexibility activation costs and minimization of active power losses.

2. The proposed DMS tool is extended for planning purposes such as to propose efficient sizing and placement of BESS solutions (i.e., distributed or centralized).

3. An analytical study is conducted to compare the alternatives among OLTC, BESS, active network management, or their coordinated operation for scenarios with increased DER integration.

4. A sensitivity analysis for coordinated operation between BESS and EVs exploring variable base pricing for the BESS investment and the variable price of EV flexibility.

3. Formulation of Coordinated Active Network Management Tool

This section details the statement of the proposed multi-objective unbalanced Multi-period AC-OPF (MACOPF). The formulation provides a flexible DMS framework for LV unbalanced networks. The main focus is to minimize the DSO/grid’s operation cost; hence, it is essential to minimize DER flexibility activation costs (e.g., minimize the engagement of DER on the operation), the overuse of any DSO assets (e.g., BESS, OLTC), minimize the grid losses as well as to minimize the energy costs by means of the energy imported by the upstream grid. Based on the strategy applied for each type of DER/asset, the optimization strategy (i.e., dispatchableμG, definition of EV flexibility, V2G, etc.) may be formulated respectively.

In this study, the LV distribution network is represented as a three-phase four wire unbalanced network with a multi-earthed neutral; this fact allows the application of the Kron’s reduction [43]. More analytical information on the modeling of lines and the transformer may be found in Appendix A. Each time slot is denoted byτ∈T, where the length of setTis the horizon of the desired scheduling of operation_Hτ. Let_xg,τbe the state vector for time slotτ represented by Equation (1), containing the instant angles and voltage for each bus (j∈{1,…_Nb}) and phase (ϕ∈Φ). The vector of decision variables_uτconsists of active and reactive power for each of the controllable DER (k∈{1,…_Nc}) as shown in (2). The voltage angles displacement between adjacent nodes may be considered as constant (commonly less than 10^∘ [40]); thus, the scale of the optimization problem can be reduced significantly. However, angles are analytically defined in this work due to the need of assessing phase balancing constraints:

_xg,τ=_ΘVτ∀τ∈T,_xτ∈^R(2∗3_Nb),

_uτ=_PcQcτ∀τ∈T,_uτ∈^R(2∗_Nc),

_yτ=_{pchpdchyπ,chyπ,dchytripytapϵj,conτ}

All of the auxiliary variables, forτ, are contained in the_yτ; such variables are involved in the DER or OLTC modeling as well as slackness variables to relax constraints and ensure convergence for any resolution. The setsNandJstand for the nodes (each bus has three nodes, one per phase) and the branches of the grid. Let X be the vector that contains stacked the state vectors, the decision variables, and any auxiliary variables defined asX=^{[_x1,…,_xHτ,_y1,…,_yHτ]T} . As explained also in [41], the auxiliary variables are intentionally appended as last elements of vector X to allow the flexible configurations of the stated problems (i.e., eases the calculation of derivatives and the data logging of initial points). To avoid lengthy notation on the problem statement, a symbolic variable_Xτ=(_xg,τ,_uτ,_yτ) is defined. The MACOPF is stated in Equation (4a):

minu∑τ=1_Htw1·_Πτ+_w2· _PLτΔτ+_w3· _UOLTCτ+_w4·_Tapτ+_Φp,τ,

subjected to

_Gj(_Xτ)=0,∀j,τ∈N,T,

_HSub(_Xτ)≤_Srated−_εsub,∀τ∈T,

VU_Fj(_Xτ)≤VUF¯−_εVUF,∀i,τ∈N,T,

_Vmin−_εVj,τ ≤_vj≤_Vmax+_εVj,τ ,j,τ∈N,T,

_hξ(_Xτ)=0,∀ξ,τ∈U,T,

_gξ(_Xτ)≤0,∀ξ,τ∈U,T,

where the analytical expression of the objective is the following:

_Πτ=∑k_Nb^[_ck(τ)]T·_uk,τO1:flexibility activation costs _PLτ=∑ϕ∈Φ∑i∈B^Δ_Vφ,ij2 _Rφ,ijO2:Apparent power losses _UOLTCτ=∑ϕ∈Φ^{_Vϕ,ps(τ)−_Vϕ,ps(τ−1)2}_O3:Penalize fast transitions of primary winding voltage_Tapτ=∑ϕ∈Φ^{_tcϕ (τ)−_tcϕ (τ−1)2}_O4:Cost of tap operations

and vector_ck(τ)assigns a price for the utilization of controllable DER or asset k atτ in €/kWh or €/kVArh. Multiple pricing schemes may be defined, enabling demand–response schemes. The exact form of nonlinear power flow equations is encapsulated with the nonlinear equality constraints in (4b); inequality constraints (4c) is posed to ensure that the MV/LV transformer is not loaded more than the nominal, or may provide a power cap for the LV grid energy management; (4d) inequality constraints stand for the phase balancing requirements; the boxed constraints in (4e) to maintain all nodal voltage within the preset limits. The slackness variablesεin (4c)–(4e) are applied to slightly relax the constraints and reassure the convergence of the optimizer even when the available active measures cannot strictly provide a solution into the feasible space (thus, it is enlarged). The last equality and inequality constraints (4f)–(4g) describe a generalized form pertaining to operational constraints of all controllable assets and DER.

In the objective function in Equation (4a), the term_Φp,τassigns to some of the auxiliary variables a penalty cost. Such penalty costs may be for the relaxation parametersε, as well as some penalties to prohibit the concurrence of charging and discharging as explained more analytically in the BESS model.

The mathematical form of the objective function is a combination of linear (_FL) and quadratic (_FQ ) cost functions. All the quadratic terms in the objective function are encapsulated by the following Equation (5):

F(X)=∑τ∈Tf(_Xτ)=_FL+_FQ=^cT·X+12^XT·H·X,

whereX∈^R_NX ,_NX=_Hτ·(_nx+_ny) contains all variables of the stated problem, including the auxiliary variables. Variable X can be defined in the following steps, by primarily considering an additional variable u as proposed in [44]. This u variable can be formed by applying a linear transformation_NTand a shiftr^to the extended set of the optimization variables:

r=_NT·X,

u=r−r^,

To enable flexible extension of such costs able to handle scaled linear and quadratic costs as in Equation (4a), each element of the optimizer full set of variables X (let_xi) is input as:

_xi=_{mi fa}(_ui+_zi),_xi<−_zi0,−_zi≤_xi≤_{zimi fa}(_ui+_zi),_xi>_zi

and

_fa=α,if_di=1^α2,if_di=2

where_ziprovides the option to shift the cost function,_miscales the variable_xiaccordingly, and, according to the specified input_di, the cost function may be shaped as linear or quadratic. This formulation is used to structure all the cost functions. The quadratic terms refer to the loss minimization as well as the operational cost functions for OLTC. For instance, the cost formulation: BESS, EV,μ G, is a piece-wise linear(PWL) function that is incorporated in the optimization framework through cost constrained variables. Further details regarding this mathematical formulation are given in [41].

3.1. Interior-Point Algorithm

The control framework in (4a) evidently presents a large scale nonlinear programming, since its size is dependent on the length of decision variables, which in turn are increased either with finer-grained time resolutions or with longer horizons of scheduling, i.e.,X∈^R_NX :_NX=(2∗3_Nb+2∗_Nc+2∗_NBESS+2∗_NEV+_NOLTC)·_Hτ. Additionally, the power flows are regarded as nonlinear equality constraints, enumerated as_Nnonlin=2∗3∗_Nb∗_Hτ.

The stated MACOPF is addressed using the Interior-Point (IP) primal-dual algorithm. Assume the compact formulation for the stated MACOPF by the set of Equation (10), where the set of variables is denoted byx:

minxf(x),

subject to

_gE(x)=0,

_hI(x)≤0,

_xmin≤x≤_xmax,

The corresponding Lagrangian function is given by Equation (11):

_Lp(x,λ,σ,s):=_fp(x)−^λT _gE(x)−^σT(_hI(x)+s)

where vectorsλ,σare the Lagrange multipliers for the corresponding equality (_gE(x)) and inequality constraints_hIwhich can be regarded also as equality constraints by the addition of slack variables s, such that_hI(x)−s=0. Thereafter, the augmented objective function for the unconstrained problem (penalty function)_fp(x) is defined by Equation (12):

_Lp(x,λ,σ,s):=_fp(x)−^λT _gE(x)−^σT(_hI(x)+s)⇔_Lp(x,λ,σ,s)=f(x)−^μ(k)∑j=1_Nxℓn(_xj−_xj,min)−^μ(k)∑j=1_Nxℓn(_xj,max−_xj)−^μ(k)∑j=1_Nineqℓn(_sj)−^λT _gE(x)−^σT(_hI(x)+s),

where^μ(k)stands for the logarithmic barrier parameter per iteration k. The latter is forced to monotonically reduce to 0 as iteration progresses by the minimizer.

The exact formulation of the non-convex nonlinear power flow equality constraints inflicts the certification of the second-order KKT conditions regarding local optimality of one point-solution^p* . An analytical description for the first and second order KKT conditions may be found in [45]. As per [41], particular computational techniques are proposed to remedy the singularities of the Jacobian matrix (caused by the inter-temporal dependencies of DER) which is necessary to be assessed along with the Hessian one, for the iterative process of the IP algorithm. Based on the MACOPF structure, efficient explicit calculations (exploiting sparsities) of the Jacobian and Hessian matrices are input in the optimizer not only to accelerate the convergence, but also to avoid faulty solution due to the singular Jacobian. A database is used to procure initial points (_X0) to the optimizer from past resolutions of the problem, as well as to acquire historical data regarding load and weather forecasts if they are not currently available.

3.2. Multi-Objective MACOPF Treated with a Weighted Sum Method

As presented in the problem statement, MACOPF is clearly a multi-objective optimization problem. The discipline of Pareto needs to be introduced hereby, which is posed to ensure that none of the objectives can be further improved in the search space without any major impact on the objective function. There exist several methods to support the decision maker, but, in practice, the most commonly used one is the weighted sum method; this recommends the scaling (i.e.,_ψi) and the multiplication of all the objectives with a weighting factor (i.e.,_wi) as follows:

F(X)=∑ik_wi·_Oi(X)1_ψi,suchthat_∑ikwi=1,w≥0,w=^{[_w1,…,_wk]T}_ψi=max(_Oi)

The scaling of each objective term is critical to balance their impact by balancing their order of magnitude on the aggregated objective function but is often disregarded leading to mistaken and overestimated efficient points.

3.3. Grid Constraints

3.3.1. Power Flows

The nonlinear power flows equations, at time instantτ , are formed by Equation (14), expressed as a function of phasor nodal voltage, injection from loads and the DER injection in complex form; essentially imposing that the mismatch between nodal injections and the injection from loads and DER is zero:

G(_Xτ)=_Snodal(_Xτ)+_Sload(_Xτ)−_Cp·_SDER(_Xτ)=0_Snodal(_Xτ)=[_Vτ]_{Ybus Vτ*}

where the_Sload∈^C3_Nbvector contains the complex loads for all buses of the system;_SDER∈^C_Nc the DER injections or consumptions. The sparse matrix_Cp∈^N3_Nb×_Ncis defined to map the DER net injections to the3_Nbnodes. Any(i,j)element of_Cpis zero, whereas it is one if generator j is located in bus i.

3.3.2. Voltage Unbalances

It is important to account for voltage unbalances particularly in LV grids, which inherently present unbalance nature. The integration of single phase DER may lead to much higher unbalances. The Voltage Unbalance Factor (VUF) has several definitions; however, EN 50160 standards make use of the sequence components as in Equation (15), [18]:

_VUFj[%]=|_υ2,j||_υ1,j|·100%≈|_υ2,j|·100%,

where_υ1,j,_υ2,j are positive and negative sequence components, respectively. Obviously, the technical constraints for unbalances present a non-convex nature. Nonetheless, the phase balancing constraint can be easily convexified by the accurate approximation that the magnitude of positive sequence components are closely to 1 p.u. as per [46]. For a node j, the phase balancing constraint is given by Equation (16):

_VUFj,τ−_εVUFj ≤2%∀j∈N,τ∈T,

_εVUFj ≤0is an auxiliary variable relaxing the balancing constraint and ensure convergence of MACOPF.

3.3.3. Line Congestion Management

The connection of several EVs along the distribution feeder may increase the peak load profile leading to line congestions. The technical constraint to manage the branch currents is applied at each phase. Let_Ib(τ)∈^C3·3·_Nbrepresent the line currents forτ slot. Exploiting the definition of the Bus-Injection to Branch-Current (BIBC) matrix as proposed in [47], the phase currents at each branch can be derived from Equation (17):

_Ibτ =BIBC·_Iτ⇔=BIBC·_{Ybus Vτ}⇔=_Ym·_Vτ,

where_Ym represents the modified admittance matrix that maps the nodal current injection vector to the respective branch current injection. The group of nonlinear inequality constraints provided in (18) are posed to ensure line congestion management. The analytical contribution of these constraints into the Jacobian and Hessian is input to the optimizer and is indicatively formulated in Appendix B:

_Ij,τ−_εL≤_Ij¯,∀j∈J,τ∈T.

3.4. OLTC and DER Operational Model

3.4.1. On-Load Tap Changer (OLTC) Model

The OLTC mechanism is considered to be connected at the primary winding of the MV/LV transformer. The primary side is connected through a branch line to the slack bus through a fictitious line impedance (i.e., Z-Thevenin) representing this way the upstream connected MV distribution network, as represented in Figure 1. The incorporation of the OLTC introduces discrete decision variables to determine the tap-positioning at each time step of the horizon. This would resort the MACOPF’s formulation to a Mixed-Integer Nonlinear (MINLP), which is generally classified as an NP-hard problem. The MINLP—especially non-convex problems—are characterized by the challenge of handling the nonlinearities in addition to the combinatorial nature posed by integer decision variables [48]. The continuous relaxation of non-convex MINLP is itself a global optimization problem, thus likely to be NP-hard [49]. In the literature, some approaches have been proposed to treat the discrete nature of tap positions by introducing continuous decision variables as in [38,39,50]. Nevertheless, none of these works couple the stages among them and subsequent time slots (i.e., multi-period optimization) or provide the option to follow any technical limitations as a maximum number of tap changes.

A three-stage resolution is hereby proposed to avoid the introduction of integer variables. In the first stage, the tap changer decision variables are treated as a continuous set[_tpa ̲,_tpa ¯] . Those decision variables are tracked with heuristic variables that follow the tangent between two-subsequent time slots. The algorithmic diagram of the proposed scheme for the OLTC is illustrated in Figure 2. The corresponding mathematical expressions that connect the tap-positioning decision variables with primary and secondary winding voltage are in Equation (19):

_VLV=_Vps−Δs·_tpϕ ,∀τ,ϕ∈T,Φ

where_VLV,_Vpsare the voltage magnitudes on secondary and primary winding accordingly andΔsis the % resolution of tap-step_tpϕ . The constraints have to be applied to all phases since many OLTC provide the option to have different tap-position per phase, yet their control is mechanically coupled to shift them all simultaneously.

The post-processing routine on the second stage lies on the round-off rule described in Algorithm 1. According to this algorithm, the continuous tap-position set is projected to the closest integer variable, whereas a tangent rule is used to enhance the algorithm and avoid excessive tap changes.

Algorithm 1: OLTC round-offξ —rule based on [50]; in this study,ξ=0.5 .

Sustainability 12 03332 i001

An optional scheme is included within the second stage relying on a minimizer to adapt the decision taken regarding the tap-position in case of particular technical limitations regarding the daily number of tap changes. This optimizer is set in Equations (20) and (21):

min_ud∑τ^{_ud(τ)−ta_pτ2}

such that

∑τ∥_ud(τ)−_ud(τ−1)∥≤Δtap¯.

The optimizer targets to reduce the distance between vectors_ud(τ)andta_pτ, such that the maximum number of tap-positioning posed for the horizon of the optimization—Δtap¯ . This formulation clearly presents a quadratic function with non-convex constraints. The problem is resolved by using a state-of-the art optimizer called NOMADS [51]. This solver implements a Mesh Adaptive Direct Search algorithm (MADS), which is capable of dealing with non-smooth objective functions and constraints since it resorts to black-box optimization, avoiding the evaluation of costly derivatives.

The post-processing of second stage outputs the discrete tap-positions, which compose inputs for the third stage. The final operational decisions are decided in the third stage where any additional control actions may be determined. The vector^X*contains control actions derived from the first stage extracting the OLTC decision variables, and this is then used as an initial point for the last MACOPF resolution.

The characteristics of the OLTC equipment used for this case study are given in Table 1. This equipment may be installed to retrofit an existing transformer, adding the OLTC capability.

3.4.2. Battery Energy Storage System (BESS)

The BESS model is a first order model, where two distinct auxiliary variables participate in the BESS state equations and operational constraints—one for the discharging_pdch≥0,_pdch∈_R+, while the charging mode_pch≤0,_pch∈_R−. Any losses occur in each mode of operation are associated with charging and discharging efficiencies (_ηch,_ηdch)._E0is defined the initial (i.e.,τ= 0) stored energy at the BESS. The remaining stored energy of a BESS at one time stepτ can be calculated by Equation (22), which clearly bundles the instant energy state with the former one:

E(τ)=E(τ−1)−Δτ_ηch1_ηdchp(τ),wherep(τ)=_pch(τ)_pdch(τ).

Within the proposed optimization framework and the subsequent participation of BESSs in the power flow equations (i.e., in_SDER), a primary decision variable per BESS defines the scalar variable for active power injections_PBESS:

_PBESS(τ)=_pch(τ)+_pch(τ),

_pch¯≤_pch(τ)≤0,

0≤_pdch(τ)≤_pdch¯,

SoC̲≤SoC(τ)≤SoC¯,

SOC(τ)=E(τ)_Erated

E(0)=E(_Hτ).

The constraints (23a)–(23e) are settled∀τ∈T. The constraints (23b)–(23d) limit the maximum charging and discharging power as well as the minimum and maximum State-of-Charge (SoC)—defined in Equation (23e)—according to the BESS’s technology and characteristics. The last constraint (23f) imposes that BESS’s ending energy state should be equal to the initial stored energy; thus, the BESS does not get fully discharged.

To avoid simultaneous charging and discharging of the BESS, a penalty cost is assigned with each auxiliary decision variables_pch,_pdch, both of which should be greater—at least one order—than the cost of use of BESS(_cBESS)itself, i.e.,_PBESS.

Based on the European Commission’s study in [26], where several scenarios for Li-ion BESS costs are concerned depending on different market growth indexes, Table 2 presents the selected characteristics for BESS used in this study. A base price is selected for the year 2025, assuming a moderate adoption of Li-ion BESS by the market. The Levelized Cost of Energy (LCOE) for BESS is also calculated to assign it with the operational costs.

3.4.3. Electric Vehicles

The same mathematical formulations as the BESS’s are settled to simulate the EV operation and technical constraints, with the difference that, in the energy state Equation (22), a term_Aτ=_ytrip(τ)·_Etr(τ)is deducted whenever there is a trip occurrence. The variable_{[ytrip]ntr×Hτ}captures the temporal occurrence of a trip combined with the energy consumed_Etralong that period. The variable_ntrsimply refers to the expected number of trips per EV. The combination of those variables creates a mapping of flexibilities (i.e., to charge or discharge based on EVs’ availability, when not used). The Vehicle-to-Grid (V2G) may be considered as an additional mode of operation for the EVs, described by the discharging decision variable.

One can define the energy state of an EV at time periodτthrough_EEV(·)∈Rthat is temporally coupled with the prior period’s energy state and the decision of charging/discharging set-point. For one EV_j , the energy state equation (22) can be recasted in a matrix format (Equation (24)) towards evolution of time as a linear combination of the initial stored energy_EEVj=^{[_EEVj(0),⋯,_EEVj(_Hτ)]T}, representing all energy states:

_EEVj=I⋮I_E0j+Λ0⋮⋱Λ⋯Λ_pEVj(1)⋮_pEVj(_Hτ)−_ytripj·_Etrj,

whereΛ=[diag{_ndch}diag{1/_nch}]·Δτ.

The EVs are modeled to emulate realistic behavior, using a mobility routine as explained analytically in [41]. Further assumptions and details about the EVs are given in Section 4, particularly in regard to the consideration of mobility model and habit of trips. The flexible use of smart charging operation is illustrated in Figure 3b, where some charging energy slots are shifted at later hours. In Figure 3a, the cost function definition is illustrated, where the dumb charging is added on the objective function as a negative cost (i.e., profit), whereas the decision to use smart charging decreases proportionally this profit. The cost function is deployed as a piecewise linear function with a CCV. Note that, for periods when no dumb charging does occur, a V (not purely symmetric since use of V2G is considered more expensive) cost function is considered with its vertex at (0,0).

3.5. Optimal Sizing and Placement of BESS

The optimal sizing and placement of the BESS is hereby derived using an extended formulation of MACOPF as presented in (25):

minu,_zjw1·_O1+_w2·_O2+_wplan·∑j∈D_zj

subject to:

k·_zj≤_pch(τ)≤0,0≤_pdch(τ)≤k·_zj,η̲_Ej(τ)_zj≤SoC(τ)≤η¯_Ej(τ)_zj.

Variable_zjis added referring to the sizing of the BESS to be installed and the setDcontains all candidate nodes for the investment of a single BESS owned by the DSO. Accordingly, the cost associated with the investment is added on the objective function. The c-rating of the BESS is denoted with k, and the under investment BESS is considered to be a three-phase system. The parametersη̲,η¯define the technical parameters for the minimum and maximum SoC for the technology of the BESS to be invested.

4. Case Study Synopsis

The validation of the proposed coordinated control of LV operation takes place for an IEEE LV benchmark network [52]. This LV grid—in Figure 4—presents the same technical characteristics with the benchmark with the difference that the MV/LV transformer 250 kVA, 20/0.4 kV, since only one feeder is regarded.

The MACOPF formulation is used to obtain the most suitable placement and sizing of the BESS for the examined grid. The results obtained, given the requirement of an energy BESS application—c-Rating 0.5—suggest the installation of a 90 kWh (round-trip efficiency considered to be 0.8) at node 460, and alternatively for distributed BESS solution one additional placed at node 666. The optimal sizing and placement problem performed co-optimizing the planning and operational services of the BESS (i.e., without the coordination with other assets or DER) for the 8-day dataset in 1-hour resolution (i.e., 384 time slots) for the mixed DER integration scenario 02; the same sizing and placement obtained are also used for all scenarios. In this particular planning stage, the BESS solution is used as the sole option to resolve the issues.

The definition of the examined scenarios is in Table 3. The PV and EV integration refer as a percentage proportional to the number of units installed out of the total 55 consumers. Note that all DER are connected to the grid as single-phase units connected in the same phase as the respective end-user. Analytical information regarding the point of connections and the characteristics of DER can be found in [41]:

A data-set has been created for this study case to emulate and induct the results for yearly analysis. A data-pool of yearly load profiles (found in repository [53]) is statistically processed by using a k-means algorithm. Several sets of load profiles are aggregated and normalized to represent MV/LV substation’s profiles. This data are then properly clustered into seasons, weekends and weekdays. Regarding the clustering each centroid is considered as the component-wise median of that cluster, let it_ci . Each datapoint is accordingly clustered in Equation (27):

d(x,c)=∑i=1_nt|_ci−_xi|.

The representative data-set is thereafter composed of eight days, two per season (one representing weekdays and another one for weekends). Each of those is selected based on the centroid metric derived by k-means. The normalized aggregated at the substation level profiles are illustrated in Figure 5. To reproduce the load profiles, the inversed cumulative Gaussian distribution function (^Φ−1) is set with maximum standard deviationσ=0.08and median the value of the centroid at each datapoint. Therefore, each point of the load profile i at instantτ is calculated from Equation (28):

_Pi(τ)=^Φ−1(_ci(τ),σ)·_Pi,rated.

The Gaussian copula method is used to generate N temporal scenarios PV solar profiles, encapsulating the seasonal dependence as proposed in [54].

Prior to the presentation of the analytical case study, an assessment of incremental DER integration (EVs and PV) on operational grid constraints is conducted. In this part of the study, no controls are deployed; though any voltage issues, voltage unbalances, and any line congestions are recorded along the yearly data-set. For the following analysis, it should be noted that, for each scenario of DER integration, five days of the same season and type of day (i.e., weekday or weekend) are considered, and the resulting metrics (i.e., voltage magnitudes and unbalances) are averaged. The collected information is exported by sequential three-phase unbalanced power flow for averaged 30-minutes profiles. Obviously, the impact of DER in higher resolutions (e.g., order of minutes) may be more intensive (i.e., leading to more severe technical problems) particularly for phase unbalances and voltage issues. All the subsequent Figure 6, Figure 7 and Figure 8 represent the per season impact of integrating residential PV and EV. For all the incremental scenarios, random values of PV units (1.7, 2.7, 3.7 kW) and EVs (charging power outlet of 3.7 and 7.4 kW) are assigned accordingly.

The evolution of maximum and minimum voltage magnitudes met in the grid for integration of PV and EV appear in Figure 6a,b, respectively. One can notice the impact of inversed power flows due to microgeneration which lifts up both maximum and minimum voltage magnitudes. During summer periods, higher overvoltage are faced; meanwhile, higher overvoltage are more common during weekdays, since the loading conditions are lighter during sunny hours (see Figure 9a). In Figure 6a, the upper x-axis describes the PV integration correlated with the peak load as a percentage. For the considered load profiles, significant voltage increase (up to 1.07 p.u.) effects close to 50% of PV integration. In higher time resolution, voltage may experience, instantaneously, values higher than 1.1 p.u.; however, this study examines only a 30-minute average profile to determine maximum and minimum voltage limits as well to identify other technical bottlenecks on the safe operation of the grid. Regarding the EV integration, the minimum voltage up to 0.915 p.u. voltage appear with a number 40 EVs along the end-users. During weekends, due to lesser EV mileage, the demand for charging appears to be more limited compared to the weekdays. Furthermore, the simultaneity of charging at late evening hours (where daily peak load appears) leads to significant voltage drops.

The impact of connecting single phase DER on VUF is presented in Figure 7a,b. In all cases where no DER is installed, there are quite balanced conditions only up to 1.3%. The installation of microgeneration leads to an increase of VUF up to 2.1% during summer periods, which is justified due to lighter load conditions and higher PV generation. The connection of EVs increases a high amount of imbalance particularly in winter periods, where the loading conditions are intensive reaching up to 2.2%. The increasing integration of EV obviously leads to more severe phase unbalancing conditions. Being in accordance with IEC standards necessitates that VUF has to be less than 2% for 95% of the week [55]; hence, for this presented analysis—that is performed for 30-minute time slots—the VUF threshold is considered at 1.8%.

The connection of multiple PV units does not arise any line congestions up to 33% of their integration. Higher PV integration, though, results in excessive reversed power along the upstream grid, particularly overloading line 280–566 up to 0.98 p.u. of line currents. The line loading conditions are much more severe concerning the EV integration as illustrated in Figure 8b, where line congestions may reach up to 1.5 p.u. particularly in winter periods.

4.1. Results

This section demonstrates the techno-economical results for the examined cases of Table 3 performed with the proposed MACOPF scheme, among several modes of operation as defined in Table 4. Note that the conventional mode of operation_m0, applies for controls with no particular intelligence. Such controls in this study are derived by executing sequential three-phase power flows to identify technical bottlenecks (i.e., line congestions or voltage issues). Based on an iterative process for every node with overvoltage, 5% of the produced energy is curtailed until the issue is resolved. Accordingly, EVs’ charging process is shed at the overloaded branches iteratively. The mode of operation_m1refers to the sole coordination of all types of DER orchestrated by the MACOPF. It should be mentioned that, for this section, all EVs are considered available to provide grid support (i.e., through smart-charging), once they are parked at the house premises. Modes of operation_m2and_m3extend the aforementioned coordination with the utilization of one and two BESS—owned by the DSO—installed at nodes 460 and 666, accordingly. The last mode of operation stands for the coordinated operation of the OLTC with available DER for grid support.

The assigned costs for the utilization of each of the DER’s flexibility for this analysis are presented in Table 5. Note that the active power curtailment is considered about three times larger than the price of selling the energy produced by residentialμ G in Portugal. The main concern of the proposed scheme is to maximize the integration of microgeneration by exploiting other sources of flexibility, at the stage of the scheduling of operation. Concurrently, closer to the time of the delivery, APC may be used by typical droop control functions, ensuring safe grid operation. The cost of energy not supplied is set according to reports in [56]. The respective costs for the OLTC and BESS investment as well as the corresponding cost (summarized in Table 1 and Table 2) of their utilization are discussed in the previous sections together with their models. Concerning the OLTC investment, it is considered that the transformer’s remaining lifetime—and subsequently the OLTC investment lifetime—is 15 years.

The yearly operational costs for all modes of operation are illustrated in Figure 10a, while a breakdown of seasonal cost analysis appears in Figure 10b. One can notice that, in most scenarios,_m0leads to higher operational costs than any coordinated operational scheme, apart from case 01 (i.e., where all modes are comparable). The increased connection of EVs in scenarios 02 and 04 leads to very high operational costs for_m0, due to the need for EV shedding in order to avoid branch congestions.

Commenting on case 01, where solely PV integration is regarded, one can notice that the cheapest operational modes of operation appear to be either the coordinated curtailment of PVs or the investment of OLTC. The curtailment ofμGshould be followed by a compelling compensation fee as applied hereby; otherwise, such schemes should adopt fairness strategies. In this particular case, MACOPF on_m1provides the optimal dispatch of active and reactive power leading to lower costs than in_m0 . Indicatively, the curtailed energy for case 01 is in Figure 11a. Both cases_m2−_m3lead to higher operational costs due to the topological distance of the BESS from the most problematic nodes (i.e., with overvoltage). Additionally, the BESS constraint for cyclic charging substantially increases the operational cost of BESS’s usage, since the absorbed energy—during sunny periods—will have to be consumed in other time slots, whether it is needed or not. For instance, BESS can be incorporated in MACOPF for extended energy management applications, where BESS are allowed to participate in the electricity market. In that case, grid operational costs for BESS usage would be more encouraging results, since the overall coordination of grid operation would follow the market price accordingly. More distributed BESS solutions may treat overvoltage byμGmore efficiently, depending on their spatial distribution along the grid. Nonetheless, the economical results between_m3–_m4are closely equivalent since the tool decides the utilization of BESS placed at node 460 in both cases. Clearly, the OLTC presents very low annual operation costs in the vicinity of 350–700 € (i.e., depending on the remaining lifetime of the retrofitted transformer) for case 01, where solely PV are installed.

The installation of EVs in scenario 03 (no PVs) is followed by undervoltage and line congestions as analyzed previously. The yearly estimated operational costs due to the shedding of loads (_m0 ) reach up to 4.3 k€. The coordinated smart charging presents a much cheaper alternative of 1.6 k€. This price may slightly fluctuate considering the uncertainty of EV users that are willing to charge under this regime. It should be reminded that smart-charging hereby is strictly regarded when EVs are not in trip progress and refers to any deviation from the expected–dumb charging profile. An example from a typical winter day used in the simulations is in Figure 11b. The dumb charging profile for case 03 results in overloaded lines up to 1.09 p.u., a fact that is addressed with the coordinated smart charging in_m1 . It can be observed that some charging profiles are shifted in slots with lesser loading conditions for the grid, which is the early morning hours as it appears in Figure 11b. In the same figure, it can be noticed that there is no V2G participation since the high EV availability concerned is capable of addressing the technical issues. In this scenario, the BESS coordinated operation presents much lower grid operational costs close to 750€ and about 700€ for_m2and_m3 , accordingly. The OLTC in this scenario—considering the 15 years lifetime of the transformer—presents annual costs close to 1.9 k€, due to the need for coordination with smart charging which pertains to 80% of this cost. Therefore, it is observed that the OLTC can only reduce the lines’ overload by lowering its taps; thus, the transformer’s secondary winding voltage is lifted up reducing the line currents analogously. This branch current reduction in this case of 5% (i.e., 1.09 p.u. reduced to 1.037 p.u.) on the most congested branch. From Figure 10b, it can be observed that the highest share of the operational costs comes from winter and autumn period, when the used data-set appears to have the peak demands.

Both scenarios 02 and 04 examine the integration of mixed DERs, considering extensive integration of PV and EV units. Case 02 is a PV rich scenario with the 30 single phaseμG, which corresponds to 49% of the peak demand of the grid. The last scenario 04 refers to a higher EV integration (35 EVs) and a 33% of the peak demand installed PV units. In both cases, the_m0 leads to high operational costs due to the need of EV shedding as well as an amount of 1.1 k€—case 02—for APC. Indicatively, Figure 12 presents analytically the control decisions derived from MACOPF for each mode of operation. It can be observed that along_m1 , smart charging schedules shift some of the EVs during sunny period hours (i.e., 10.00–14.00), while some other EVs are further charged in the beginning of the day before the trip occurrences. The V2G mode of operation takes place not only to reduce line currents, but also to create available charging slots to be used during sunny periods with expected overvoltage. All EVs are constrained to keep their SoC at the end of the day equal to the one at the beginning of the simulation. Therefore, EVs that are parked at home all day present high availability, which is observed to be used in this way (i.e., charge during periods with high solar irradiance and discharge to avoid line congestions). Regarding OLTC operation, as illustrated in Figure 12d, it ends during the end of each day at tap –2, which is due to the loss minimization term. The OLTC acts in such way to increase the secondary winding voltage in order to avoid any voltage drops and minimize the active power losses as well. Obviously, the addition of some more EVs (case 04) results in a very abrupt increase of operational costs for_m0, a fact that is connected to the extensive line congestions that occur.

The proposed OLTC equipment is allowed to be set offline in different tap positions per phase to treat unbalances. The performed analysis considers 30-minute data resolution, a fact that may underestimate the real-time conditions regarding phase unbalances. Therefore, the OLTC may need further engagement with flexible DER to deal with phase unbalances; hence, higher operational costs may be foreseen in such a case.

4.1.1. Minimization of Active Power Losses

The annual active power losses without DER integration are estimated to be 1.5% to 2.5%, calculated through sequential power flow executions. From Figure 13, it can be observed that there is an increase in line losses due to the extensive DER integration, particularly when EVs are considered. The charging of the EVs does most likely occur during peak load period in the afternoon, increasing notably the loading along the distribution lines and to some extent the line losses (see Figure 13 reaching close to 6%).

In all cases and the subsequent modes of operation towards the resolution of technical bottlenecks (i.e., voltage magnitudes, voltage unbalances, and line congestions) line losses are also reduced compared to the scenario where no controls are applied. Particularly, modes_m1–_m4based on the proposed MACOPF control framework do further reduce the line losses due to the involved objective term. The coordinated operation of the OLTC with DER appears to perform the most efficient measure in the direction of losses minimization. Additionally, experimenting (i.e., by assigning higher values to_w2) with the weighted terms among the objective terms in modes_m1–_m3does not impact significantly the control decisions and the losses. The latter can be justified due to the fact that further minimization of active power losses, substantially, precedes more engagement of DER flexibility. Nonetheless, involving more flexibility in the scheduling of operation for the the losses minimization cannot be justified economically.

The exploitation of the OLTC brings greater benefit to the objective, particularly for the minimization of active power losses, since the low cost per tap operation allows the optimizer to derive taps-down, lifting the primary side to higher voltage.On the other hand, the combination of BESS or smart charging cannot justify their excessive usage towards further minimization of losses at the same order of reduction compared to the OLTC.

4.1.2. Sensitivity Analysis on BESS and Smart-Charging Coordination

A sensitivity analysis is presented hereby to observe the evolution of the estimated DSO annual operational costs for different prices of smart charging and different LCOE for the BESS. This section provides a comprehensive comparative analysis not only in the reflection of the pricing of EV flexibility—i.e., EV available for smart charging—and the BESS, into the DSO equivalent annual cost of operation, but also signifies the importance of enabling EV smart charging towards the reduction of operational costs. Furthermore, one scenario concerns the V2G mode of operation in the sensitivity analysis. Along the evolution of BESS and EV pricing and their availability (i.e., spatio-temporal), the resulting coordination is recorded through the proposed control scheme.

In each case, the range of BESS prices lies within [200–380] €/kWh (i.e., corresponding LCOE [0.05, 0.0925] €/kWh) according to the study in [26], for base year 2025. For the purpose of this study, scenario 02 is considered regarding the integration of DER. The cost of utilizing EV flexibility is considered in the range [0.05, 0.25] € for each kWh of shifted charging slot.

On the data-points in Figure 14, the coordination of EV and BESS is displayed as a share of the annual DSO operation costs. As BESS price increases, the EVs’ flexibility is further used following a quasi-linear dependence and vice versa. Note that, particularly in Figure 14a, where only 20% of the connected EVs are considered as flexible, even when BESS cost is higher than the EV flexibility, the limited EV availability forces the resolution of technical bottlenecks mainly with the use of BESS. On the other hand, one can notice that, in Figure 14b, the higher availability of EVs—60%–(which implies better spatio-temporal distribution) leads to reduced DSO costs in most cases (see data-points in Figure 14b). The fact that EVs are distributed along the grid is foreseen to be very efficient to address any technical bottlenecks arising from PV and EV integration. The latter can be further observed when all EVs are considered as flexible in Figure 14c, where the effect on the cost reduction is more intensive since the plane follows—an affine—curvature. The last sensitivity analysis considers V2G in constant pricing at 0.35 €/kWh, combined with 100% EV availability for smart charging. The fact that V2G is considered in the framework essentially allows the utilization of EVs that do not proceed with a trip (i.e., and then they have to get charged, since they are constrained to recharge in the simulation day at least the energy used). The effective coordination of the BESS, EV smart charging, and V2G leads to notably lesser grid operational costs compared to other cases. Note that the maximum cost observed when BESS and EV are assigned with the highest price values is 1.6 k€, when the resulting maximum cost in the other planes reaches close to 2.5 k€.

5. Conclusions

This work presents a control framework that derives the coordinated management for multiple DSO assets and DER, in a way to ensure safe grid operation pertaining admissible voltage magnitudes, phase balancing, and avoiding line congestions. The overall scheme is formulated as a nonlinear multi-objective program resolved with a primal-dual interior point algorithm. A three-stage technique is proposed to incorporate the OLTC in the decision-making process. The control framework is advanced and adapted to propose efficient decisions for the placement and sizing of BESS.

The study shows that OLTC appears to be the most efficient option to treat overvoltage when high PV integration is encountered, considering loss minimization. Nonetheless, phase unbalances may occur that could be treated by coordinating with other DER, or the installation of BESS. It was concluded that, for long-term phase imbalances, the OLTC can be also setup in offline mode to proper tap positions per winding—hence reducing the unbalances and the need of coordination with other assets or DER. The extensive integration of EVs (more than 20 EVs in the examined case) cannot be accommodated only by the optimal operation of the OLTC. In the presented study the OLTC is capable of reducing the overloading to 25–32% of the most congested branches; however, coordination with EVs is deemed necessary to respect all technical constraints. Depending on the selected OLTC technology and whether the transformer can be retrofitted—as examined in this study—the OLTC is foreseen to be most efficient to address reversed power flows effects and the subsequent overvoltage. On the other hand, BESS’s solution is very dependent on the expected costs along the evolution in the next decade. Considering a moderate BESS cost (i.e., regarding the expected adoption of BESS in the market) presents comparable results—or better results in some cases—with OLTC. In mixed DER scenarios (i.e., PV and EV), BESS coordinated with DER outperforms the compared modes of operation, presenting the lowest DSO annual equivalent operational costs.

A sensitivity analysis demonstrated the DSO’s annual costs of operation considering different costs for BESS and EV flexibility utilization. The main outcome of this sensitivity study shows the importance, due to their spatio-temporal distribution, of the active participation of DER (i.e., mainly for EVs with V2G hereby) in the grid operation and the significant reflection on cost reduction.

Investment cost (€)	cinv	7.000
Step Voltage (%)	Δs	(up to 3) hereby constant at 2
Min/Max tap-position	tap̲/tap¯	(up to±9)±2
Min/Max voltage (p.u.)	Vps̲/Vps¯	1.1/0.9
Maintanance-free operations	Noltc	700.000
Approximated Cost per Tap (€)	ctap	0.01

	2025 (reference year)
Price (euro/kWh) *(includes costs of investment)	290
Cycles DoD at 80% in lifetime	5000
LCOE calculation (€/kWh)	0.0725

	Case 01	Case 02	Case 03	Case 04
PV [nr of PV units]	30	30	0	20
EV [nr of EVs]	0	30	30	35

Operational Mode	Conventional Operation –No Smart Controls Applied– (m0)	DER Optimal Coordination (m1)	BESS Coordinated with DER (m2)	Distributed BESS Coordinated with DER (m3 )	Coordination of OLTC with DER
OLTC	✕	✕	✕	✓	✓
BESS	✕	✕	✓	✕	✓
Smart Charging	✕	✓	✓	✓	✓
Vehicle to Grid	✕	✓	✓	✓	✓
μ G Active Power Curtailment	✓	✓	✓	✓	✓
μ G Reactive Power Dispatch	✕	✓	✓	✕	✕
Load Shedding	✓	✕	✕	✕	✕

	B	Cost (€/kWrh-kVArh)
Cost of Active Power Curtailment	capc	0.30
Cost Smart Charging	capc	0.15
Cost of V2G	cV2G	0.35
Cost of Energy Not Supplied Lines	cENS	3

Author Contributions

Conceptualization, K.K., I.M., J.L.D.-G., and N.H.; Formal analysis, K.K.; Investigation, K.K.; Methodology, K.K.; Software, K.K.; Supervision, I.M., J.L.D.-G., H.L., N.S., and N.H.; Validation, K.K. and J.L.D.-G.; Writing-original draft, K.K.; Writing-review and editing, I.M., J.L.D.-G., H.L., N.S., and N.H. All authors have read and agreed to the published version of the manuscript.

Funding

This project has received funding from the European Unions Horizon 2020 research and innovation program under the Marie Skodowska-Curie Grant No. 675318 (INCITE).

Acknowledgments

The authors would like to thank Luís Marques for his significant support on the statistical processing of the data-sets.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Distribution Network Models (Lines and Transformer)

The representation of line model as well as the the interconnection of the LV with upstream distribution grid are in Figure 1. All nodes of the grid in this study have three terminals, each of which represents the phasesa,b,c. The voltage magnitude for node j is given by the real vector_vj=^{[_vj,a,_vj,b,_vj,c]T}, whereΦ={a,b,c}is the set of available phases in the distribution grid. Accordingly, the voltage angles by the real vector_ϑj∈^R3.

A connection between buses j and k is mathematically represented by a square symmetric matrix_zk,m∈^{C_Φk,m×_Φk,m}(e.g., Kron's reduction, the analytical form would explicitly include the neutral and the earth conductor), where_Φk,mis the number of phases of interconnected nodes k and m. The active conductors (i.e., the three-phases, neutral follows the reduction) present coupling amongst them; hence, the[_zk,m]has off-diagonal elements different from 0, as well as the corresponding self-inductances. The admittance matrix (_Ybus∈^C3_Nb×3_Nb ) defines the topological structure and the connectivity among nodes of the distribution network. The line shunt admittances for the distribution lines in LV grid can be neglected [57]). Consequently, the element_{Ykpk ,mpm} of_Ybuswhich refers to the connection between phase_pkof bus k and phase_pmof m can be expressed as:

(_{Ykpk ,mpm} )=-1_(zk,m)pk,pmifk≠m,_∑n1_zk,nifk=m.

The ideal voltage source is used to impose constant phase in the slack bus_θsource=^{0,2π3,-2π3T}as well as to represent a stiff busbar according to its equivalent Thevenin impedance (i.e., calculated by the short-circuit power) placed in series. The_Vsourcefor power flow applications does represent the slack bus; instead, within the proposed optimization, it follows the voltage derived by OLTC (_Vps).

The distribution transformer (MV/LV) is also included in the system's admittance matrix according to its type (i.e., mainly the type of connection of the windings). The contribution of the transformer can be represented and included in the_Ybus of the network by constant impedances (i.e., for steady state analysis) representing . Based on the transformer's winding connection, the admittance matrices for the distribution transformer can be found in [58].

Appendix B. Calculation of Derivatives

As this constraint presents a nonlinear one, the analytical first and second derivative will be given to the optimizer. For the first order derivative of the branch current injections, the derivation with regard to the vector X is given in (A3) and the subsequent derivative for a particular set of state and decision variable of time-step_xτk is given in (A4). More analytically, it is:

∂_Ib,τk∂_Θτk =_Ym·∂_Vτk ∂_Θτk =j_Ym·[_Vτk ]

∂_Ib,τk∂_Vτk =_Ym·∂_Vτk ∂_Vτk =_Ym·[_Vτk ]^{[_Vτk ]-1}

Obviously, the partial derivatives of branch currents for period_τk(_Ibτk) with regard to_Pg,τkand_Qg,τkare zero entries to the Jacobian matrix, together with all the rest of the derivations that correspond to decision and control variables from other time-steps (e.g.,∂_Ibτk/∂_x1=0).

Accordingly, the second-order derivatives for the vector that correspond to the branch currents for time-step_τk will be structured in proportion to the Lagrangian multiplier for inequality constraints, by the following matrix (A5):

∂_Ib,τk∂X=∂_Ib,τk∂_x1⋯∂_Ib,τk∂_xτk ⋯∂_Ib,τk∂_xHt ∂_Ib,τk∂_xHt ∂_Ib,τk∂_y1⋯∂_Ib,τk∂_yHt

∂_Ib(_τk)∂_xτ=∂_Ib(_τk)∂_Θτ∂_Ib(_τk)∂_Vτ∂_Ib(_τk)∂P_gτ∂_Ib(_τk)∂Q_gτ

σ∂2 Ibτk∂xτk2=σ∂∂xτk Ibτk∂xτk T=Ib,τk Θτk Θτk Ib,τk Θτk Vτk 00Ib,τk Vτk Θτk Ib,τk Vτk Vτk 0000000000

where

_{Ib,τk Θτk Θτk} (σ)=∂∂_Θτk j[_Vτk ]_YmTσ=-[_YmTσ][_Vτk ]

_{Ib,τk Vτk Θτk} (σ)=∂∂_Θτk [_Vτk ]^{[_Vτk ]-1} _YmTσ[V]=j[_YmTσ][_Vτk ]^{[_Vτk ]-1}=-j_{Ib,τk Θτk Θτk} σ^{[_Vτk ]-1}

_{Ib,τk Θτk Vτk} (σ)=∂∂_Vτk ^{∂_Ib,τk∂_Vτk T}σ=j[_YmTσ][_Vτk ]^{[_Vτk ]-1}=_{Ib,τk Vτk Θτk} (σ)

_{Ib,τk Vτk Vτk} (σ)=∂∂_Vτk ^{∂_Ib,τk∂_Vτk T}σ=0