A Horizon Optimization Control Framework for the Coordinated Operation of Multiple Distributed Energy Resources in Low Voltage Distribution Networks

Abstract

In recent years, the installation of residential Distributed Energy Resources (DER) that produce (mainly rooftop photovoltaics usually bundled with battery system) or consume (electric heat pumps, controllable loads, electric vehicles) electric power is continuously increasing in Low Voltage (LV) distribution networks. Several technical challenges may arise through the massive integration of DER, which have to be addressed by the distribution grid operator. However, DER can provide certain degree of flexibility to the operation of distribution grids, which is generally performed with temporal shifting of energy to be consumed or injected. This work advances a horizon optimization control framework which aims to efficiently schedule the LV network’s operation in day-ahead scale coordinating multiple DER. The main objectives of the proposed control is to ensure secure LV grid operation in the sense of admissible voltage bounds and rated loading conditions for the secondary transformer. The proposed methodology leans on a multi-period three-phase Optimal Power Flow (OPF) addressed as a nonlinear optimization problem. The resulting horizon control scheme is validated within an LV distribution network through multiple case scenarios with high microgeneration and electric vehicle integration providing admissible voltage limits and avoiding unnecessary active power curtailments.

1. Introduction

Nowadays, an increasing number of small-scale units, typically referred to Distributed Energy Resources (DER), is connected along the Low Voltage (LV) distribution networks posing several technical challenges, whilst bringing novel and diverse opportunities. Most commonly there is already a large share of injected power at the distribution level by the renewable energy merely based on solar energy through Photovoltaic (PV). The connection of such resources at the LV grid and end-users’ premises is foreseen to increase substantially in the close future with small rooftop installation usually coupled with Battery Storage System (BESS), controllable loads (e.g., Electric Water Heaters) and Electric Vehicles (EV). Therefore, it shall be critical to exploit the DER controllability and active participation through demand response schemes in order to support or even improve techno-economically the operation of the distribution networks operation [1,2]. For instance, DERs might be used to provide technical support to tackle voltage or congestion problems delivering profits to the end-consumers accordingly.

Traditionally, LV networks used to be the most passive circuits within the power systems, since power flows were solely headed from distribution transformers to consumers without the operation of automation elements [3]. In particular, the entire segment from the secondary substation and its downstream connected LV networks is very often not monitored nor controlled [4].

The Distribution System Operators (DSOs) address such technical challenges by increasing the observability and controllability of the grids, envisioning the active management of the DERs for ancillary services, throughout new operation stages [5,6]. Such attributes of advanced control and monitoring techniques do typically refer to Advanced Distribution Management Systems (A-DMS) [7]. The active management of distribution networks through the engagement of DERs in the operation of the grid is regarded to occur with the provision of flexibilities services such as active and reactive power control (i.e., inverter based control). The smart grid deployment, as an alternative paradigm for the operation of distribution networks, envisions the active management of DER taking advantage of advanced control infrastructures and communications through demand side management schemes. Advanced control methodologies need to be implemented to determine control actions related to controllable DER, which can techno-economically improve distribution networks’ operation delivering benefits to residential users.

Virtual Power Plants and aggregation of flexibilities in higher levels for market participation (i.e., in balancing market) and provision of ancillary services through DER has been foreseen in [8,9]. In several European countries already, flexible loads are incentivized in order to defer network investments or decrease congestions by load shifting [1].

A particular concern in recent research works regards the potential flexibility provided by the DER connected along the LV distribution level to address voltage problems. In past years, research was focused on active power curtailment following local droop based rules or even combined with reactive control [10,11,12]. Moreover, self-consumption is a common practice adopted by the DSOs lately, in order to avoid voltage rise during the peak period of PV generation. In several European countries (e.g., Belgium, Denmark and the Netherlands), residential PV self-consumption measures based on net metering schemes target matching the endogenously generated power with local consumption [13]. On the contrary, Portugal and Germany promote lower remuneration for energy produced by microgeneration, thus attracting instantaneous consumption. Towards the path to maximizing renewable generation into distribution networks, the focus in research remains in controlling the microgeneration itself. A distributed scheme with more sophisticated rules to mitigate overvoltages due to a large integration of PVs was proposed in [14]. In addition, real-world LV four-wire distribution networks are in practice fairly unbalanced, since single-phase grid elements (e.g., end-users, micro-generation and EVs) do impact the voltage, not only of the connected phase, but also of the other two phases due to the neutral-point shifting phenomenon [15]. Consequently, local droop based controls via single-phase PVs might be insufficient for voltage regulation in unbalanced grids; hence, the deployment of centralized optimization schemes can tackle such issues since topological (i.e., making use of the most efficient controllable DER or asset) and phase coupling considerations among phases can be regarded [11,14].

Interest has also been attended for the efficient integration and exploitation of distributed BESS [16,17]. Recently, BESS has been introduced by electric utilities to accommodate the increased generated power by solar energy in LV grids, though the deployment of residential BESSs has been limited up to the previous years, due to the relatively high capital cost of such devices. Lately, the continuous reducing cost of batteries in addition to the rising electricity costs and incentives for investments in storage [18]. According to Directive 2009/72/EC [19], the utilization of energy storage systems by grid operators is very limited at present; meanwhile, unbundling requirements for DSOs under EU directives do not allow energy storage units to be directly owned, or controlled by them. Concurrently, the growing number of BESS owned by residential consumers is likely to undermine the current business model of the electric utilities [20]. The following trend aims to maximize the revenue brought to the consumer in particular when home energy management systems are utilized to optimize the local generation and consumption.

Other research works have proposed advanced controlling more DER types such as EV and Controllable Loads (CL) for the mitigation of overvoltages or line congestions by [21,22,23,24,25,26]. A centralized control scheme for the voltage regulation and the mitigation of high unbalance instances is proposed [27], the efficiency of which is compared with typical local control droops. An extension of the same authors provide a framework for the coordination of an On Load Tap Changer (OLTC), installed at the secondary substation, with BESS and controllable microgeneration [20].

Optimal Power Flow (OPF) is widely applied as a tool within DMS application for the planning and operation of the power systems. Clearly, OPF problems are deemed challenging since they require solving of non-convex problems. Nonconvexity stems from the nonlinear relationship between voltages and the complex powers consumed or injected at each node [28]. Further adaptations and assertions have to regarded for power flow equations in particular for LV grids as they present purely unbalanced loading conditions and mainly resistive line characteristics. The widely used DC power flow methodologies in transmission grid studies, but cannot be applied due to the higher R/X ratios [29]. The application of non-convex and nonlinear AC power flows in an OPF framework possibly leads to computational complexity according to [30]. Therefore, in literature, convex relaxations are settled, based on e.g., semidefinite relaxations [31]; such approaches explore solutions that are globally optimal for the original problem in many practical cases, leading though non-valid solutions in some cases [28]. Further computational complexity can be certainly added to the OPF formulation if it is settled for multiple periods leading to multi-period OPF.

Recent works have dealt with proposing efficient linearizations to resort tractable multi-period OPF [6,17]. For instance, authors in [17] take advantage of the linearization to reduce the convergence of the programming and utilize it for planning of the distribution network. In [6], a step forward advances the multiperiod OPF framework incorporating uncertainties brought by forecasts through chance constrained optimization. Nevertheless, both works do not address the three-phase nature and the subsequent unbalances may arise in LV distribution networks. In this work, a three-phase multi-period OPF based on the exact (i.e., nonlinear) AC power flows is proposed, incorporating multiple DER within the operation of the distribution grid.

The main contributions of this paper can be summarized as follows:

• A decision tool which provides support to the DSO for the minimization of the operational costs based on the coordinated operation of multiple DER. The tool is capable of mitigating the regulating of nodal voltages, minimizing curtailments of active power by the microgeneration, ensuring nominal rated power for the secondary for all time instances. Multiple active measures are posed based on different DER technologies.
• A three-phase multi-period OPF framework based on the exact formulation of the AC power flow equations. The overall problem is resolved through a nonlinear optimization problem addressed interior-point method where efficient explicit calculation for the gradients of the constraints and the Hessian of the Lagrangian are proposed leaning on sparsities.
• Analytical inter-temporal constraints (i.e., providing the limitations of each type of DER) and the counterpart inter-temporal cost dependencies are discussed with their subsequent burdens. In particular, a technique is proposed to address singularity of Jacobian matrix (i.e., of the nonlinear problem) induced by the inter-temporal constraints.

This paper is structured as follows: Section 2 describes the mathematical models of the power system as well as the temporal models deployed for the DER. Section 3 discusses the proposed coordinated control scheme, which is mathematically stated through a three-phase Multiperiod AC-OPF (MACOPF). An analytical discussion regarding the resolution of the nonlinear programming through the Interior-Point (IP) primal-dual algorithm takes places in the same section. Section 4 and Section 5 present the study for the validation of the proposed control scheme. The final section provides the conclusions.

2. Proposed Models for Distributed Energy Resources and Distribution Network Models

2.1. Distribution Network Models (Lines and Transformer)

The LV distribution network is assumed to be comprised by a typical three-phase four wire unbalanced network with a multi-earthed neutral; this fact allows the application of the Kron’s reduction [32]. All buses are considered to have three terminals, where each one represents the phase connection point: $a,b,c$. The voltage magnitude for bus j is given by the real vector $v_j\in {\mathbb{R}}^3$, $v_j={[v_{j,a},v_{j,b},v_{j,c}]}^T$, and, accordingly, the voltage angles by the real vector ${\vartheta }_j\in {\mathbb{R}}^3$.

The connection between buses k and m is a square (e.g., Kron’s reduction) symmetric matrix $z_{k,m}\in {\mathbb{C}}^{{\mathrm{\Phi }}_{k,m}\times {\mathrm{\Phi }}_{k,m}}$, where ${\mathrm{\Phi }}_{k,m}$ the number of phases interconnected nodes k and m. The active conductors (e.g., the three-phases, neutral follows the reduction) present coupling, amongst them the $\left[z_{k,m}\right]$ has off-diagonal elements different from 0. The admittance matrix ($Y_{bus}\in {\mathbb{C}}^{3N_b\times 3N_b}$) determines the overall topological structure of the distribution network, where assuming no shunt admittances modeled for the distribution lines (i.e., negligible according to [33]), the element $Y_{k_{p_k},m_{p_m}}$ of $Y_{bus}$ which refers to the connection between phase $p_k$ of bus k and phase $p_m$ of m can be expressed as:

$$(Y_{k_{p_k},m_{p_m}})=\left\{\begin{array}{cc}-\frac{1}{{\left(z_{k,m}\right)}_{p_k,p_m}}& \mathrm{if}k\ne m,\\ {\sum }_n\frac{1}{z_{k,n}}& \mathrm{if}k=m.\end{array}\right.$$

The distribution transformer (MV/LV) according to its type most commonly delta–wye grounded, wye–delta, wye–wye, open-wye–open-delta, delta–delta can be represented and included in the $Y_{bus}$ of the network by constant impedances (i.e., for steady state analysis). For each type of transformer configuration, the admittance matrices for the distribution transformer can be found in [34]. Building the admittance matrix of the transformer is commonly known that it is not always invertible other than wye-g-wye-g; hence, an addition of a fictitious small admittance from the isolated transformer sides to the ground remedies the issue [35].

The interconnection of the LV grid with the MV grid is represented with an ideal voltage source that is assigned with the voltage vector $V_{source}={[1,1,1]}^T$, which is also considered as the reference bus with angles $\phase{V_{source}}={\left[0,\frac{2\pi }{3},-\frac{2\pi }{3}\right]}^T$. The ideal voltage source is connected in series with proper impedances that can be assessed according to the MV grid’s short circuit power (Figure 1).

2.2. Battery Energy Storage System (BESS)

The BESS model is based on a first order battery model. Two distinct auxiliary variables are settled as power injections. The positive term refers to the discharging mode of operation $p_{dch}\ge 0$, $p_{dch}\in {\mathbb{R}}_{+}$, while the charging of the storage unit is negative $p_{ch}\le 0$, $p_{ch}\in R_{-}$. This model captures the losses during charging and discharging modes, through the corresponding efficiencies (${\eta }_{ch},{\eta }_{dch}$). ${\mathcal{E}}_0$ is considered the initial (i.e., $\tau$ = 0) stored energy of the BESS. The available energy capacity of a BESS at time step $\tau$ can calculated by Equation (1), which bundles the instant energy state with the former one:

$$\begin{array}{c}\mathcal{E}\left(\tau \right)=\mathcal{E}(\tau -1)-\mathrm{\Delta }\tau \left[\begin{array}{cc}{\eta }_{ch}& {\eta }_{dch}\end{array}\right]p\left(\tau \right),\end{array}$$

(1)

where

$$\begin{array}{c}p\left(\tau \right)=\left[\begin{array}{c}{p}_{ch}\left(\tau \right)\\ p_{dch}\left(\tau \right)\end{array}\right].\end{array}$$

The energy state for the BESS for the last step will be accordingly defined as a linear combination with the previous states of its energy. As it will be presented in Section 3.4, such inter-temporal couplings led to problematic conditions for the resolution of an optimal control, hence special treatment is proposed. In the proposed optimization framework, as a primary decision variable for each BESS is considered, its power injection $P_{BESS}$, which should be subjected to the equality constraint (2a) for each instant $\tau$, followed by some operational constraints for both operation modes as follows:

$$\begin{array}{c}{P}_{BESS}\left(\tau \right)=p_{ch}\left(\tau \right)+p_{ch}\left(\tau \right),\end{array}$$

(2a)

$$\begin{array}{c}\overline{p_{ch}}\le p_{ch}\left(\tau \right)\le 0,\end{array}$$

(2b)

$$\begin{array}{c}0\le p_{dch}\left(\tau \right)\le \overline{p_{dch}},\end{array}$$

(2c)

$$\begin{array}{c}\underline{SoC}\le SoC\left(\tau \right)\le \overline{SoC},\end{array}$$

(2d)

$$\begin{array}{c}SOC\left(\tau \right)=\frac{\mathcal{E}\left(\tau \right)}{{\mathcal{E}}_{rated}}.\end{array}$$

(2e)

The constraints (2a)–(2e) are posed $\forall \tau \in \mathcal{T},\mathcal{T}:=\{1,\dots ,H_{\tau }\}$ and $H_{\tau }$ the last instant defining the horizon of the optimization. The constraints (2a)–(2e) pose the technical constraints for the BESS charging and discharging power. Accordingly, the State-of-Charge (SoC)—defined in Equation (2e)—is constrained based on the BESS’s characteristics. To avert simultaneous charging and discharging of the BESS, a penalty cost is assigned with each auxiliary decision variables $p_{ch},p_{dch}$, both of which should greater—at least one order—than the use of the BESS $\left(c_{BESS}\right)$ itself, i.e., $P_{BESS}$.

As an additional note, for any three-phase BESS the different phases are hereby considered to follow the same mode of operation; therefore, the mathematical expression is comprised of three single-batteries installed in the different phases, coupled with the equality constraints for their active and reactive power injections as in Equation (3):

$$\begin{array}{ccccc}{P}_{BESS}^A\left(\tau \right)& =& P_{BESS}^B\left(\tau \right)& =& P_{BESS}^C\left(\tau \right)\\ Q_{BESS}^A\left(\tau \right)& =& Q_{BESS}^B\left(\tau \right)& =& Q_{BESS}^C\left(\tau \right)\end{array}\forall \tau \in \mathcal{T}.$$

(3)

2.3. Electric Vehicles (EVs)

The EVs are structured following the same rationale as the BESS model. In this study, the EVs are considered as flexible DER according, certainly, to their availability each time. Their provided flexibility is essentially regarded to be the intervals when they are parked to their owner’s house premises. Being in this state (i.e., parked), if there is need for charge, this will be decided by the proposed control following the coordinated smart charging scheme. When the owner of an EV desires to provide a signal of flexibility, the time interval when the estimated trip will occur together with the estimated consumed energy should be dispatched to the DSO. These two signals are captured for each controllable (i.e., willing to be charged in concordance to the smart charging scheme) with $\left[y_{trip}\right]$ that is added to discharge the EV and $E_{tr}$, where $n_{tr}$ corresponds to the number of trips for an EV. The fictitious variable ${\left[y_{trip}\right]}_{n_{tr}\times H_{\tau }}$ is added to discharge the EVs during their trips.

The Vehicle-to-Grid (V2G), where the EV injects power to the grid, operation is also incorporated within the proposed EV model. Whenever the V2G mode of operation is not deemed to be followed, the $p_{dch}$ is simply constrained to zero.

One can define the energy state for each instant for one EV given by the vector ${\mathcal{E}}_{EV}\in {\mathbb{R}}^{H_{\tau }}$ recasting Equation (1), which infers to a linear combination of preceding instances inherent to the controllability that its flexibility allows, and the initial stored energy ${\mathcal{E}}_0$. The energy storage for one EV at instant $\tau$ can be calculated by Equation (4):

$$\mathcal{E}\left(\tau \right)={\mathcal{E}}_0+\underset{\mathrm{\Lambda }}{\underbrace{\mathrm{\Delta }\tau \left[\mathrm{diag}\left\{n_{dch}\right\}\mathrm{diag}\{1/n_{ch}\}\right]}}·\underset{p_{\mathrm{EV}}\left(\tau \right)}{\underbrace{\left[\begin{array}{c}{p}_{dch}\\ p_{ch}\end{array}\right]}}-y_{\mathrm{trip}}\left(\tau \right)·E_{\mathrm{tr}}\left(\tau \right).$$

(4)

For one EV, let it j, connected along the distribution network, the energy state function (4) can be rewritten in a compact matrix format capturing both operating modes where the energy stored to each EV towards the time evolution $H_{\tau }$ can expressed as the vector ${\mathbf{E}}_{\mathrm{EV}}^j={[E_{\mathbf{EV}}^j\left(0\right),\dots ,E_{\mathrm{EV}}^j\left(H_{\tau }\right)]}^T$:

$${\mathbf{E}}_{\mathrm{EV}}^j=\left[\begin{array}{c}I\\ ⋮\\ I\end{array}\right]{\mathcal{E}}_0^j+\left[\begin{array}{cc}\mathrm{\Lambda }& \mathbf{0}\\ ⋮& \ddots \\ \mathrm{\Lambda }& \cdots & \mathrm{\Lambda }\end{array}\right]\left[\begin{array}{c}{p}_{\mathrm{EV}}^j\left(1\right)\\ ⋮\\ p_{\mathrm{EV}}^j\left(H_{\tau }\right)\end{array}\right]-{\mathbf{y}}_{\mathrm{trip}}^j·{\mathbf{E}}_{\mathrm{tr}}^j.$$

(5)

2.4. Microgeneration ($\mu$G)

The microgeneration in this work is considered to be single phase inverter based installations. In case the DSO desires to incorporate the control of microgeneration in the control, the type of control pertaining the active power through Active Power Curtailment $\left(APC\right)$ and/or Reactive Power Control $\left(QR\right)$ have to be opted for the centralized controller.

Regarding the APC, the following settings define the maximum possible curtailed power as a percentage of the instant injected power (i.e., maximum curtailment $\beta$=20%), given the following rule in Equation (6):

$$\overline{p_{pac}\left(\tau \right)}=\left\{\begin{array}{cc}\beta ·p_{inj}\left(\tau \right),\hfill & \mathbf{if}{p}_{inj}\left(\tau \right)\ge \xi ·p_{rated},\hfill \\ 0,\hfill & \mathbf{else}\mathbf{.}\hfill \end{array}\right.$$

(6)

$\xi$ (in this study $\xi$=0.5) stands for a parameter which leads to control PV with larger injected power in proportion to their installed power at the instant period. The reactive power control is defined in similar way, though allowing capacitive and inductive operation (i.e., injecting and absorbing reactive power accordingly). Nevertheless, if the option of controlling the PVs in both PAC and QR mode, to avoid the nonlinear constraint inherent to the operation of the inverter; a simpler linear constraint is posed to ensure that the microgenerator’s inverter does not exceed its rated power:

$$\begin{array}{c}\overline{Q_{QR}\left(\tau \right)}=\left(p_{rated}-p_{inj}\left(\tau \right)\right)·tan\left({\theta }_{\min }\right)\\ \underline{Q_{QR}\left(\tau \right)}=-\overline{Q_{QR}\left(\tau \right)},\end{array}$$

(7)

where ${\theta }_{\min }$ is given by the minimum power factor (PF${}_{\min }$) applied; ${\theta }_{\min }={cos}^{-1}\left({\mathrm{PF}}_{\min }\right)$.

2.5. Three-Phase Power Flow

A three-phase Power Flow (PF) is implemented following the main notions described in [36]. The PF tool is incorporated in the overall proposed scheme, as an algorithmic step, either for the calculation of the initial point of the optimization, or to validate of the set-points yielded by the control scheme.

Th three-phase power is based on a Backward-Forward Sweep (BFS) technique, where, in the Backward stage, the branch current calculations occur, whilst, in the Forward Sweep stage, the nodal voltage is calculated. This method, unlikely for classical power flow methods, copes with a branch-oriented technique rather than nodal relations.

Note that this power flow algorithm presents quick convergence, i.e., iterations do not exceed 4 for tolerance convergence ${ϵ}_{\tau }=1e-4$. The performance can be further accelerated by the valid assertion that the angle displacement in LV distribution networks between adjacent nodes is fairly small [17], i.e., $\mathrm{\Delta }\theta$ leads to zero which leads to the conception in Equation (8) for the voltage drop:

$$\mathrm{\Delta }{V}_{abc}^{(\kappa +1)}=Re\{Z_{\ell }·J_{abc}^{\left(\kappa \right)}\},$$

(8)

where $Z_{\ell }$ is the corresponding impedance among the connected branches and $J_{abc}^{\left(\kappa \right)}$ is the vector for the line section currents at iteration $\kappa$. Regarding the power flow framework, it is hereby structured in such a way that each load might have a different load model among constant PQ and a constant I or constant Z model. Accordingly, the injection current at node j is given by Equation (9):

$$I_{abc}^j={(S_{abc}·dia{g}^{-1}\left(V_{L-L}\right))}^{\ast }{∥\frac{V_j}{V}∥}^{\kappa },$$

(9)

where $S_{abc}$ stands for the apparent power consumed at node j, $V_{L-L}$ the line-to-line voltage. The operator $diag(.)$ is settled as an operator that returns a diagonal vector and $\kappa$ is considered the load model parameter, which is 0 for constant PQ load, 1 for constant current and 2 in case of constant impedance.

3. Multi-Period AC-OPF (MACOPF) Formulation-Coordinated Control Scheme

The centralized coordinated management of the DER is discussed in this section through a multi-period three-phase AC-Optimal Power Flow (MACOPF) where the different periods $\tau \in \mathcal{T}$ are coupled with inter-temporal costs and the DER are assigned with inter-temporal technical constraints accordingly.

The MACOPF is stated for a horizon of operational planning $H_t$. All subsequent time steps belong to the set $\mathcal{T}:=\{1,\dots ,H_{\tau }\}$. The main objective (i.e., ${\mathcal{O}}_1$) of the scheme is to minimize the operating costs assigned with all the controllable assets providing their coordination given their availability. Penalty costs assigned to auxiliary variables described with the term ${\mathrm{\Phi }}_p$. Such penalty costs refer to relaxation of voltage bounds to ensure feasibility and thus convergence as well as penalties to avert simultaneous charging and discharging or even auxiliary variables described within Section 2.

Assume that the state vector ($x_{g,\tau }$) at the time instant $\tau$ is given by (10) and the set of decision variables $\tau$ corresponds to the vector $u_{\tau }$ comprised of active and reactive power of each controllable DER as shown in Equation (11) as well as auxiliary variables $y_{\tau }$ in Equation (12). The voltage angles can be omitted to reduce the scale of the optimization problem, since the angle displacement between adjacent nodes in LV grids is typically less than 10${}^{°}$ [17]. Nonetheless, angles are included for completeness:

$$\begin{array}{cc}\hfill x_{g,\tau }& ={\left[\begin{array}{c}\Theta \\ \mathcal{V}\end{array}\right]}_{\tau },\forall \tau \in \mathcal{T},x_{\tau }\in {\mathbb{R}}^{(2\ast 3N_b)},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill u_{\tau }& ={\left[\begin{array}{c}{P}_c\\ Q_c\end{array}\right]}_{\tau },\forall \tau \in \mathcal{T},u_{\tau }\in {\mathbb{R}}^{(2\ast n_c)},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill y_{\tau }& ={\left[\begin{array}{c}{p}_{ch}\\ p_{dch}\\ y_{\pi ,ch}\\ y_{\pi ,dch}\\ {\epsilon }_V\\ {\epsilon }_{sub}\\ y_{trip}\end{array}\right]}_{\tau },\forall \tau \in \mathcal{T},y_{\tau }\in {\mathbb{R}}^{{\ell }_y},\hfill \end{array}$$

(12)

where $N_b$ refers to the number of buses and $N_c$ the number controllable units and ${\ell }_y=(4\ast N_{BESS})+3\ast N_b+N_{ev}+3$ with $n_{BESS}$ the number of BESS and $n_{EV}$ the number of EVs. The real vector $\mathcal{V}={[v_1,v_2,\dots ,v_{N_b}]}_{\tau }^T$ corresponds to the voltage magnitudes for each bus (each bus has three terminals) at each time instant $\tau$, and respectively $\Theta$ to the voltage magnitudes. The sets $\mathcal{N},\mathcal{J},\mathcal{T}$, denote the buses ($N_b$), branches and the horizon of the multi-period scheme. Let for the optimization problem the state vector and the decision variables correspond to the respective matrices $X={[x_1,\dots ,x_{H\tau }]}^T$ and $\mathcal{U}={[u_1,\dots ,u_{H\tau }]}^T$, essentially defined as stacked vectors of each subsequent time period. For the sake of comprehension, Figure 2 presents the described structure of the optimization variables. The fact that the auxiliary variables ($y_{\tau }$) are appended as last elements of vector X eases the extension of the stated problem. Additional objective terms might be assigned in the current formulation unless there is no dependence or conflicting interest upon the aforementioned objective.

The overall control scheme can be described by the set of Equations (13):

$$\underset{u}{min}\sum _{\tau =1}^{H_t}\left\{\underset{{\mathcal{O}}_1}{\underbrace{\sum _k^{N_b}\left({\left[c_k\left(\tau \right)\right]}^T·u_{k,\tau }\right)}}\right\}\mathrm{\Delta }\tau +{\mathrm{\Phi }}_p,$$

(13)

subjected to

$$\begin{array}{ccc}\hfill G_j(x_{\tau },u_{\tau },y_{\tau })=0& \hfill & \forall j,\tau \in \mathcal{N},\mathcal{T},\hfill \end{array}$$

(13a)

$$\begin{array}{ccc}\hfill H_{Sub}(x_{\tau },u_{\tau },y_{\tau })+{\epsilon }_{sub}\le 0& \hfill & \forall i,\tau \in \mathcal{J},\mathcal{T},\hfill \end{array}$$

(13b)

$$\begin{array}{ccc}\hfill V_{\min }-{ϵ}_{\mathcal{V}}\le v_j\left(x_{\tau }\right)\le V_{\max }+{ϵ}_{\mathcal{V}}& \hfill & \forall i,\tau \in \mathcal{N},\mathcal{T},\hfill \end{array}$$

(13c)

$$\begin{array}{ccc}\hfill h_{\xi }(x_{\tau },u_{\tau },y_{\tau })=0& \hfill & \forall \xi ,\tau \in \mathcal{U},\mathcal{T},\hfill \end{array}$$

(13d)

$$\begin{array}{ccc}\hfill g_{\xi }(x_{\tau },u_{\tau },y_{\tau })\le 0& \hfill & \forall \xi ,\tau \in \mathcal{U},\mathcal{T},\hfill \end{array}$$

(13e)

where vector ${\left[c_k\left(\tau \right)\right]}^T$ includes the price of each controllable unit per time step $\tau$ in €/kWh or €/kVArh. The constraints in (13a) set the nonlinear power flow equation at each bus of the network; inequality constraint (13b) poses the technical constraint for the MV/LV transformer; the boxed constraints in (13c) to respect all nodal voltages to range strictly within the admissible bounds. The additional positive variable ${ϵ}_{\mathcal{V}}$ is used to relax the voltage constraints and avoid infeasibility. The latter is applicable substantially when the active measures are not adequate for addressing voltage problems. Accordingly, since a transformer is capable of being operated in full load conditions or slightly higher for a certain short interval, an auxiliary variable is also applied to turn these constraints less tight. Th constraints (13d)–(13e) correspond to the operational limits of the controllable DER.

In the following subsection, the resolution of the optimal control problem is analytically discussed through the optimization techniques used to address it efficiently. The proposed control scheme evidently deals with a large number of decision variables -$X\in {\mathbb{R}}^{N_X}:N_X=(2\ast 3N_b+2\ast N_c+2\ast N_{BESS}+2\ast N_{ev}+3)·H_{\tau }$-. Accordingly, the power flow through the nonlinear equations are accounted $N_{nonlin}=2\ast 3\ast N_b$. Such large optimization problems, where the structure of the Jacobian of the nonlinear constraints present very sparse blocks, should possibly reveal particular framing of the problem re-structuring it and inherently leading to improved computational efficiency. In [37], the authors restructure the multi-period OPF to a tailored approach exploiting its particularities. Hereby, particular techniques are proposed to speedup the largeness of the nonlinear optimal control programming. The explicit calculation of the Jacobian and the Hessian is proposed, taking advantage of the sparsities, as well as slight pivotal elements to the Jacobian address any singularities met by the inter-temporal constraints. The initial point $X_0$ for the optimizer is either obtained through a local database which has stored previous occurrences or by running sequentially ($H_{\tau }$) power flows. Additionally, if load and weather forecasts fail to be obtained, a worst case scenario can be interrogated by a local historical database.

3.1. Interior-Point Primal-Dual Method for the Proposed Control Scheme

The proposed control scheme based on the MACOPF is stated with the set of Equation (13). Hereby, a discussion will follow based on an interior-point method, which follows the above form of the previous section, and will involve lengthy and complicated notation. To ease the description, the control problem is reformulated in a more compact manner (i.e., grouping the equalities and inequality constraints) as proposed in the literature [38,39,40]. Note that, with bold notation, a vector form is implied. The decision variables and state vectors are simply represented by one vector $\mathbf{x}$:

$$\underset{\mathbf{x}}{min}f\left(\mathbf{x}\right),$$

(14)

subject to

$$g_E\left(\mathbf{x}\right)=\mathbf{0},$$

(14a)

$$h_I\left(\mathbf{x}\right)\le \mathbf{0},$$

(14b)

$${\mathit{x}}_{min}\le \mathbf{x}\le {\mathit{x}}_{max},$$

(14c)

where $g_E\left(\mathbf{x}\right)$ contains any type of equality constraint (i.e., linear and nonlinear) and $h_I$ any type of inequality constraint. The inequalities constraints can be introduced as equality constraints by the addition of slack variables $s_j$, such that $h_I\left(\mathbf{x}\right)-\mathbf{s}=\mathbf{0}$. Then, a penalty function introduces a new form to the initial objective function as follows:

$$\underset{\mathbf{x}}{min}{f}_p=f\left(\mathbf{x}\right)-{\mu }^{\left(k\right)}\sum _{j=1}^{N_x}\ell n(x_j-x_{j,min})-{\mu }^{\left(k\right)}\sum _{j=1}^{N_x}\ell n(x_{j,max}-x_j)-{\mu }^{\left(k\right)}\sum _{j=1}^{N_{ineq}}\ell n\left(s_j\right),$$

(15)

subject to

$$\begin{array}{c}{g}_E\left(\mathbf{x}\right)=0,\\ h_I\left(\mathbf{x}\right)+\mathbf{s}=0,\\ \mathbf{x},\mathbf{s}\ge \mathbf{0},\end{array}$$

(16)

where ${\mu }^{\left(k\right)}$ is the the logarithmic barrier parameter for iteration k that essentially reduces monotonically to 0 as iteration progresses. The non-negativity conditions at (16) are handled by incorporating them into logarithmic barrier terms. The Lagrangian (${\mathcal{L}}_p$) of the $f_p$ can be defined as:

$$\begin{array}{c}{\mathcal{L}}_p(\mathbf{x},\mathbf{\lambda },\mathbf{\sigma },\mathbf{s}):=f_p\left(\mathbf{x}\right)-{\mathbf{\lambda }}^T{g}_E\left(\mathbf{x}\right)-{\mathbf{\sigma }}^T(h_I\left(\mathbf{x}\right)+\mathbf{s})\iff \\ \begin{array}{c}\hfill {\mathcal{L}}_p(\mathbf{x},\mathbf{\lambda },\mathbf{\sigma },\mathbf{s})=f\left(\mathbf{x}\right)-{\mu }^{\left(k\right)}\sum _{j=1}^{N_x}\ell n(x_j-x_{j,min})-{\mu }^{\left(k\right)}\sum _{j=1}^{N_x}\ell n(x_{j,max}-x_j)-\\ \hfill {\mu }^{\left(k\right)}\sum _{j=1}^{N_{ineq}}\ell n\left(s_j\right)-{\mathbf{\lambda }}^T{g}_E\left(\mathbf{x}\right)-{\mathbf{\sigma }}^T(h_I\left(\mathbf{x}\right)+\mathbf{s}),\end{array}\end{array}$$

(17)

where vectors $\mathbf{\lambda }$, $\mathbf{\sigma }$ are the Lagrange multipliers for the corresponding equality and inequality constraints. The first-order Karush–Kuhn–Tucker (KKT) condition (iff $g_E,h_I$ are first order differentiable) for a local optimum point $p^{\ast }=({\mathbf{x}}^{\ast },{\mathbf{\lambda }}^{\ast },{\mathbf{\sigma }}^{\ast },{\mathbf{s}}^{\ast })$ are the following:

$$\begin{array}{cccc}{\nabla }_x{\mathcal{L}}_p({\mathbf{x}}^{\ast },{\mathbf{\lambda }}^{\ast },{\mathbf{\sigma }}^{\ast },{\mathbf{s}}^{\ast })& =& {\nabla }_x{f}_p\left(x\right)-\mathbf{\lambda }{\nabla }_x{g}_E{\left(\mathbf{x}\right)}^T-\mathbf{\sigma }{\nabla }_x{h}_I{\left(\mathbf{x}\right)}^T-\mathcal{D}& =\mathbf{0},\\ {\nabla }_{\mu }{\mathcal{L}}_p& =& \mu {\left[\mathbf{s}\right]}^{-1}e-\lambda & =\mathbf{0},\\ {\nabla }_{\lambda }{\mathcal{L}}_p& =& g_E\left({\mathbf{x}}^{\ast }\right)& =\mathbf{0},\\ {\nabla }_{\sigma }{\mathcal{L}}_p& =& h_I\left({\mathbf{x}}^{\ast }\right)-{\mathbf{s}}^{\ast }& =\mathbf{0},\\ \mathcal{D}=\mathbf{\mu }{[\mathbf{x}-{\mathbf{x}}_{min}]}^{-1}\mathbf{e}+\mathbf{\mu }{[{\mathbf{x}}_{max}-\mathbf{x}]}^{-1}\mathbf{e},\end{array}$$

(18)

where $\mathbf{e}$ is the appropriate size vector with all entries equal to one.

The first-order KKT conditions are necessary and sufficient for global optimality for a convex problem when the Linear Constraint Qualification (LICQ) holds [41]. In the proposed control scheme, the non-convex nature of the—exact—nonlinear power flow necessitates the verification second-order KKT conditions to certify the local optimality of $p^{\ast }$. The second-order conditions can be found analytically in the literature [38], since here the first-order will be analytically discussed due to the examination of linear dependence of the inter-temporal constraints introduced by DER.

The IP algorithm at iteration k requires the solution of the following nonlinear system:

$$\left[\begin{array}{cccc}{\nabla }_x^2{\mathcal{L}}_p& {\nabla }_x{g}_E\left(\mathbf{x}\right)& {\nabla }_x{h}_I\left(\mathbf{x}\right)& 0\\ {\nabla }_x{g}_E{\left(\mathbf{x}\right)}^T& 0& 0& 0\\ {\nabla }_x{h}_I{\left(\mathbf{x}\right)}^T& 0& 0& I\\ 0& 0& I& {\nabla }_s^2{\mathcal{L}}_p\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }\mathbf{x}\\ \mathrm{\Delta }\mathbf{\lambda }\\ \mathrm{\Delta }\mathbf{\sigma }\\ \mathrm{\Delta }\mathbf{s}\end{array}\right]=-\left[\begin{array}{c}{\nabla }_x{f}_p\left(x\right)-\mathbf{\lambda }{\nabla }_x{g}_E{\left(\mathbf{x}\right)}^T-\mathbf{\sigma }{\nabla }_x{h}_I{\left(\mathbf{x}\right)}^T\\ \mu {\left[\mathbf{s}\right]}^{-1}e-\lambda \\ g_E\left(\mathbf{x}\right)\\ h_I\left(\mathbf{x}\right)-\mathbf{s}\\ \left[\mathbf{z}\right]\mathbf{x}-\mu \mathbf{e}\end{array}\right].$$

(19)

The KKT system in Equation (19) is nonlinear and its solution most commonly in the literature [42] is addressed by applying the Newton’s method. In the proposed control scheme, the gradients for the nonlinear constraints and the Hessian of the Lagrangian are explicitly calculated. In case where these derivations are not provided, approximations based on finite differences are typically applied [43,44].

3.2. Gradients of Nonlinear Constraints and Hessian of Lagrangian

The gradient (see also Appendix B) of the objective function, the Jacobian of nonlinear constraints and Hessian of the Lagrangian are implicitly provided to the optimizer, by expanding the calculations presented in [45]. On this point, it is important to be stated that the voltage vectors are expressed using polar coordinates: expressing voltage in rectangular coordinates eliminates trigonometric functions from the PF equations [46]. Nevertheless, in [42,47], a benchmarked comparison of both types of coordinates present the same order of computational performance as well as an equivalent number of iterations for convergence for a typical OPF. On one hand, rectangular coordinates can ease the process of first and second order gradients leading to quadratic and constant forms; both types cannot avoid the nonlinear equalities (and inequality) constraints and form a convex region. Concurrently, rectangular coordinates may provide slightly faster evaluation of particular gradients and Hessian, but the voltage bounds are handled as functional bounds in many OPF problems.

Therefore, the complex voltage vector might be denoted by $V\in {\mathbb{R}}^{3N_b}$. The element at bus j at phase a is ${\upsilon }_{j,a}=\left|{\upsilon }_{j,a}\right|e^{j{\theta }_{j,a}}$. The derivation with the state vectors for instant period $\tau$ followed, obviously, for the other periods are 0 entries-:

$$\begin{array}{cc}{V}_{{\mathrm{\Theta }}_{\tau }}=\frac{\partial V_{\tau }}{\partial {\mathrm{\Theta }}_{\tau }}=j\left[V_{\tau }\right]& V_{{\mathcal{V}}_{\tau }}=\frac{\partial V_{\tau }}{\partial {\mathcal{V}}_{\tau }}=\left[V_{\tau }\right]{\left[V_{\tau }\right]}^{-1}:=\left[E\right].\end{array}$$

(20)

The analytical AC power flow equations over all periods of the horizon window can be derived by the resolution of Equation (21). The operator ⊙ is used for element-wise matrix product. Complex number equations are not addressed by state-of-the art optimizers and only in a few cases yield faster solutions [48]. Thus, a segregation is proposed—as shown in Equation (22)—which introduces the power flows through the vector $G\left(X\right)$, due to their complex nature, which essentially cannot be posed at the optimization stage:

$$\begin{array}{c}\begin{array}{cc}{G}_c\left(X\right)=S_{bus}\left(V\right)+S_d-C_g{S}_g& ,G_c:{\mathbb{C}}^{n_b}\to {\mathbb{C}}^{n_b}\end{array}\\ S_{bus}^{\tau }\left(V_{\tau }\right)=V_{\tau }\odot {\left(I_{inj}^{\tau }\right)}^{\ast }=V_{\tau }\odot {(Y_{bus}·V_{\tau })}^{\ast },\end{array}$$

(21)

$$\begin{array}{cc}G\left(X\right)=\left[\begin{array}{c}\Re \left\{G_c\left(X\right)\right\}\\ \Im \left\{G_c\left(X\right)\right\}\end{array}\right]& ,G:{\mathbb{R}}^{2n_b}\to {\mathbb{R}}^{2n_b}.\end{array}$$

(22)

The current bus injection $I_{inj}$ appears in the power flow Equation (21). It would be useful for the power flow expressions derivation to present the corresponding for the current injections:

$$\begin{array}{c}\frac{\partial I_{inj}^{\tau }}{\partial x_{\tau }}=\left[\begin{array}{cccc}\frac{\partial I_{inj}^{\tau }}{\partial {\mathrm{\Theta }}_{\tau }}& \frac{\partial I_{inj}^{\tau }}{\partial {\mathcal{V}}_{\tau }}& 0& 0\end{array}\right],\\ \begin{array}{cc}\frac{\partial I_{inj}^{\tau }}{\partial {\mathrm{\Theta }}_{\tau }}=Y_{bus}\frac{\partial V_{\tau }}{\partial {\mathrm{\Theta }}_{\tau }}\stackrel{\left(20\right)}{=}j{Y}_{bus}\left[V_{\tau }\right],& \frac{\partial I_{inj}^{\tau }}{\partial {\mathcal{V}}_{\tau }}=Y_{bus}\frac{\partial V_{\tau }}{\partial {\mathcal{V}}_{\tau }}\stackrel{\left(20\right)}{=}{Y}_{bus}\left[E\right].\end{array}\end{array}$$

(23)

The first and second derivatives for the nonlinear constraints, which substantially refer to the power flow equation, will be based on the introduced vector $G_c\left(X\right)$:

$$G_X=\frac{\partial G}{\partial X}=\left[\begin{array}{c}\frac{\partial \Re \left\{G_c\left(X\right)\right\}}{\partial X}\\ \frac{\partial \Im \left\{G_c\left(X\right)\right\}}{\partial X}\end{array}\right]=\left[\begin{array}{cccccccc}\frac{\partial G}{\partial x_1}& \frac{\partial G}{\partial x_2}& \cdots & \begin{array}{|c|}\hline \frac{\partial G}{\partial x_{H_t}}\\ \hline\end{array}\cdots & \frac{\partial G}{\partial x_{H_t}}& \frac{\partial G}{\partial y_1}& \cdots & \frac{\partial G}{\partial y_{H_t}}\end{array}\right],$$

(24)

$$\frac{\partial G}{\partial x_{\tau }}=\left[\begin{array}{cccc}{G}_{{\mathrm{\Theta }}_{\tau }}& G_{v_{\tau }}& G_{P{g}_{\tau }}& G_{Q{g}_{\tau }}\end{array}\right].$$

(25)

The first order partial derivatives are presented for the defined $G_c$, which can thereafter appended in Equation (25):

$$\begin{array}{c}{G}_{c,{\mathrm{\Theta }}_{\tau }}=\frac{\partial S_{bus}\left(V_{\tau }\right)}{\partial {\mathrm{\Theta }}_{\tau }}=[I_{inj}^{\tau }]\frac{\partial V_{\tau }}{\partial {\mathrm{\Theta }}_{\tau }}+[V_{\tau }]\frac{\partial {\left(I_{inj}^{\tau }\right)}^{\ast }}{\partial {\mathrm{\Theta }}_{\tau }}\stackrel{\left(20\right)}{=}j[V_{\tau }]{\left([I_{inj}^{\tau }]-Y_{bus}[V_{\tau }]\right)}^{\ast },\end{array}$$

(26)

$$\begin{array}{c}{G}_{c,{\mathcal{V}}_{\tau }}=\frac{\partial S_{bus}\left(V_{\tau }\right)}{\partial {\mathcal{V}}_{\tau }}=[I_{inj}^{\tau }]\frac{\partial V_{\tau }}{\partial {\mathcal{V}}_{\tau }}+[V_{\tau }]\frac{\partial {\left(I_{inj}^{\tau }\right)}^{\ast }}{\partial {\mathcal{V}}_{\tau }}\stackrel{\left(20\right)}{=}[V_{\tau }]{\left([I_{inj}^{\tau }]+Y_{bus}[V_{\tau }]\right)}^{\ast }{[{\mathcal{V}}_{\tau }]}^{-1},\end{array}$$

(27)

$$\begin{array}{cc}{G}_{c,P_g^{\tau }}=-C_g,& G_{c,Q_g^{\tau }}=-j{C}_g,\end{array}$$

(28)

where ${\left(C_g\right)}_{N_b\times N_c}$ stands for the injection connectivity matrix that each element $(i,j)$ is one if at bus $i_{th}$ controllable asset $j_{th}$ is connected, else the element is zero. It is obvious that the partial derivatives of $G_c\left(x_{\tau }\right)$ with respect to other $x_{g,{\tau }_2},u_{g,{\tau }_2}$—with ${\tau }_2\ne \tau$—results in zero entries. The $G_X$ is therefore formed by two stacked matrices, which present sparse block diagonalities. The corresponding partial derivatives with respect to the auxiliary variables will be zero entries as well.

The second derivatives for the the complex power injections are necessary only for the assessment of the Hessian of the Lagrangian, which appears in iteration of the KKT-system—Equation (18). As it can be observed, the Hessian matrix of the Lagrangian can be given by Equation (1):

$${\mathcal{H}}_p(\mathbf{x},)={\nabla }_{XX}{L}_p={\nabla }_{XX}f\left(X\right)+{\nabla }_{XX}{g}_E\left(\lambda \right)+{\nabla }_{XX}H\left(X\right).$$

(29)

The second order derivatives for the complex power flows are calculated in proportion to the instant $\lambda$. The derivative is provided discretized in two parts:

$$\begin{array}{ccc}\frac{{\partial }^2{G}_c\lambda }{\partial x_{\tau }^2}& =& \frac{\partial }{\partial x_{\tau }}\left(\mathbf{\lambda }{G}_{c,x_{\tau }}^T\left(X\right)\right)\\ & =& \left[\begin{array}{cccc}{G}_{c,{\mathrm{\Theta }}_{\tau }{\mathrm{\Theta }}_{\tau }}& G_{c,{\mathrm{\Theta }}_{\tau }{\mathcal{V}}_{\tau }}& 0& 0\\ G_{c,{\mathcal{V}}_{\tau }{\mathrm{\Theta }}_{\tau }}& G_{c,{\mathcal{V}}_{\tau }{\mathcal{V}}_{\tau }}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].\end{array}$$

(30)

A brief presentation of all the subsequent expressions follows:

$$\begin{array}{ccc}{G}_{c,{\mathrm{\Theta }}_{\tau }{\mathrm{\Theta }}_{\tau }}& =& \frac{\partial }{\partial {\mathrm{\Theta }}_{\tau }}\left(j{\left([I_{inj}^{\tau }]-[V_{\tau }]Y_{bus}^T\right)}^{\ast }[V_{\tau }]{\mathbf{\lambda }}_{\tau }\right)\\ & =& \underset{A_1}{\underbrace{{[V_{\tau }]}^{\ast }\left({\left(Y_{bus}^{\ast }\right)}^T[V_{\tau }][{\mathbf{\lambda }}_{\mathbf{\tau }}-\left[{\left(Y_{bus}\right)}^{\ast }[V_{\tau }]{\mathbf{\lambda }}_{\mathbf{\tau }}\right]\right)}}+\underset{A_2}{\underbrace{{\mathcal{C}}_2\left(Y_{bus}^{\ast }{[V_{\tau }]}^{\ast }-{[I_{bus}]}^{\ast }\right)}}.\end{array}$$

(31)

Accordingly, the $G_{c,{\mathcal{V}}_{\tau }{\mathrm{\Theta }}_{\tau }}$ can be calculated as:

$$G_{c,{\mathcal{V}}_{\tau }{\mathrm{\Theta }}_{\tau }}=j\mathcal{B}(A_1-A_2)={\left(G_{c,{\mathrm{\Theta }}_{\tau }{\mathcal{V}}_{\tau }}\right)}^T,$$

(32)

$$G_{c,{\mathcal{VV}}_{\tau }}=\mathcal{B}(\mathcal{C}+{\mathcal{C}}^T)\mathcal{B},$$

(33)

where

$$\begin{array}{ccc}\mathcal{B}\hfill & =\hfill & {\left[V_{\tau }\right]}^{-1},\hfill \\ {\mathcal{C}}_2\hfill & =\hfill & \left[{\lambda }_{\tau }\right]\left[V_{\tau }\right],\hfill \\ \mathcal{C}\hfill & =\hfill & {\mathcal{C}}_2{\left(Y_{bus}\left[V_{\tau }\right]\right)}^{\ast }.\hfill \end{array}$$

The overall assessment of the second order gradient can be induced to a simple routine reducing computational effort and memory allocation by saving certain matrices presented above.

The subsequent gradients of the objective functions can be assessed for the proposed scheme since the $f\left(X\right)={\sum }_{\tau \in \mathcal{T}}f\left(x_{\tau }\right)$ substantially corresponds to a linear combination of convex functions. The first-order gradient of the objective will be comprised of constant and null entries since the costs are linear functions, as described in Section 3.5. By extension, the second-order derivatives of the objective function will be a null matrix.

If a solution of the barrier problem satisfies the primal-dual Equation (19) of the nonlinear KKT system, its solution may be approximated by an iteration of Newton’s method. The search direction can be obtained as a solution of the linearization of the KKT system, which is presented in Section 3.3.

3.3. Solution of Karush–Kuhn–Tucker Equations

Hereby, the barrier solution is presented for the optimal control problem (14), assuming that the decision variables x are positive (i.e., to avoid lengthy notation). The dual variables can be introduced as:

$$z_i:=\frac{\mu }{x_i}$$

(34)

$${\left[\begin{array}{ccccc}\mathcal{H}& \mathbf{0}& -J_E^T& -J_I^T& -I\\ \mathbf{0}& \left[\mathbf{\lambda }\right]& \mathbf{0}& \left[\mathbf{s}\right]& \mathbf{0}\\ J_E& 0& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ J_I& -I& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \left[\mathbf{z}\right]& \mathbf{0}& \mathbf{0}& \mathbf{0}& \left[\mathbf{x}\right],\end{array}\right]}^{\left(k\right)}{\left[\begin{array}{c}\mathrm{\Delta }\mathbf{x}\\ \mathrm{\Delta }\mathbf{s}\\ \mathrm{\Delta }\mathbf{\lambda }\\ \mathrm{\Delta }\mathbf{\sigma }\\ \mathrm{\Delta }\mathbf{z}\end{array}\right]}^{\left(k\right)}=-\left[\begin{array}{c}{\nabla }_x{\mathcal{L}}_p\\ {\nabla }_{\lambda }{\mathcal{L}}_p\\ {\nabla }_{\sigma }{\mathcal{L}}_p\\ {\nabla }_s{\mathcal{L}}_p\end{array}\right],$$

(35)

where $\mathcal{H}={\nabla }_{xx}{\mathcal{L}}_p$, $J_E={\nabla }_x{g}_E\left(\mathbf{x}\right)$ and $J_I={\nabla }_x{h}_I\left(\mathbf{x}\right)$. The system (35) can be further simplified—based on Gaussian elimination—by eliminating the last row. Therefore, with this elimination, it will be:

$${\left[\begin{array}{cccc}\widehat{\mathcal{H}}& \mathbf{0}& -J_E^T& J_I^T\\ \mathbf{0}& \left[\widehat{\mathbf{\lambda }}\right]& \mathbf{0}& -I\\ J_E& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ J_I& -I& \mathbf{0}& \mathbf{0}\end{array}\right]}^{\left(k\right)}{\left[\begin{array}{c}\mathrm{\Delta }\mathbf{x}\\ \mathrm{\Delta }\mathbf{s}\\ \mathrm{\Delta }\mathbf{\lambda }\\ \mathrm{\Delta }\mathbf{\sigma }\end{array}\right]}^{\left(k\right)}=-\left[\begin{array}{c}{\nabla }_x{\mathcal{L}}_p+{\left[\mathbf{x}\right]}^{-1}{g}_E\left(\mathbf{x}\right)\\ {\left[\mathbf{s}\right]}^{-1}(-\mu \mathbf{e}+\mathbf{s}\lambda )\\ -g_E\left(\mathbf{x}\right)\\ -h_I\left(\mathbf{x}\right)\end{array}\right]$$

(36)

where $\widehat{\mathcal{H}}=-{\mathcal{H}}^k-{[{\mathbf{x}}^k]}^{-1}[{\mathbf{z}}^k]$ and the diagonal matrix $[\widehat{\mathbf{\lambda }}]=-({\mathbf{s}}^k)[\mathbf{\lambda }]$. The updates for the dual variables can be obtained from Equation (37):

$$\mathrm{\Delta }{\mathbf{z}}^k={\left[{\mathbf{x}}^k\right]}^{-1}(\mu \mathbf{e}-\left[{\mathbf{z}}^k\right]\mathrm{\Delta }{\mathbf{x}}^k)-{\mathbf{z}}^k.$$

(37)

The overall resolution of the KKT system provides the search direction set for each subsequent iteration. An important condition for the Netwon step is that the Jacobian $J_E$ is a non-singular matrix. The latter implies properties which are assigned with the constraint qualification (CQ). The QP is a critical condition that needs to be assessed along with the KKT conditions. More analytically, the LICQ necessitates that the gradients of the equality constraints and any active bound constraints (i.e., binding constraints) which are linearly independent. If this does not hold, the KKT system cannot be resolved. In [41], the authors show that, for any local optimizer, the KKT conditions with LICQ satisfied can ensure the generic existence of the Lagrangian multipliers. The inter-temporal couplings among different periods $\tau$ in some cases lead to singularities for the Jacobian of the nonlinear equalities. Section 3.4 thoroughly discusses an efficient manner to avert such issues.

3.4. Intertemporal Couplings and Singular Jacobian

In the proposed multi-period OPF, the inclusion of inter-temporal constraints (e.g., mainly due to BESS and EVs) in most cases lead to the singularity of Jacobian matrix ($J_E$). Whenever the gradients of the active constraints are linearly dependent, the consequence is that the Jacobian matrix for the first-order optimality conditions will be singular. This can be observed for the presented KKT systems (36) when there are some—assuming a set of j-binding constraints (from the $h_I$ inequalities) ${\mathbf{s}}_{\mathbf{j}}=0$, while ${\mathbf{\sigma }}_{\mathbf{j}}\ne \mathbf{0}$. Therefore, if additionally there are gradients of the respective gradients $\nabla g_E,\nabla h_I$ are dependent when constraints are active the last three rows of the iteration matrix (36) will have dependent rows. An additional problematic condition might appear whether the binding conditions are linearly dependent; then, the Jacobi matrix is again singular according to [49].

The singularity issue to some extent occurs due to the structure of the control scheme, which essentially has no unique set of Lagrange multipliers corresponding to the dependent binding constraints. The particular case of BESS and the subsequent problems are analytically discussed in [49,50]. In these works, techniques are proposed to address the singularities; in [49], the authors simply suggest the elimination of the linearly dependent binding constraint on the Jacobian matrix is singular. Nonetheless, none can guarantee that the no constraint violation will take place and meanwhile it is tailored to mitigate certain singularities. The same authors propose further methods based on either a Moore–Penrose pseudoinverse of the Jacobian or by adding the standby losses.

Hereby, a simple technique will be presented which is model-free, which aims to correct the singularity introducing pivotal changes within the Jacobian matrix based on notions presented in [38]. To avoid the failure of singularity, a small pivot element (i.e., shadowed elements) might be added whenever the issue arises

$${\left[\begin{array}{cccc}\widehat{\mathcal{H}}& \mathbf{0}& -J_E^T& J_I^T\\ \mathbf{0}& [\widehat{\mathbf{\lambda }}]& \mathbf{0}& -I\\ J_E& \mathbf{0}& {\mathbf{\delta }}_{\mathbf{\pi }}·\mathit{I}& \mathbf{0}\\ J_I& -I& \mathbf{0}& \mathbf{0}\end{array}\right]}^{\left(k\right)}{\left[\begin{array}{c}\mathrm{\Delta }\mathbf{x}\\ \mathrm{\Delta }\mathbf{s}\\ \mathrm{\Delta }\mathbf{\lambda }\\ \mathrm{\Delta }\mathbf{\sigma }\end{array}\right]}^{\left(k\right)}=-\left[\begin{array}{c}{\nabla }_x{\mathcal{L}}_p+{\left[\mathbf{x}\right]}^{-1}{g}_E\left(\mathbf{x}\right)\\ {\left[\mathbf{s}\right]}^{-1}(-\mu \mathbf{e}+\mathbf{s}\lambda )\\ -g_E\left(\mathbf{x}\right)-{\mathbf{\delta }}_{\mathbf{\pi }}·\mathbf{\lambda }\\ -h_I\left(\mathbf{x}\right)\end{array}\right],$$

(38)

where ${\delta }_{\pi }$ stands for a fairly small positive number. The latter correction takes place in each iteration that the Jacobian is tracked as rank-deficient, which implies the linear dependency among constraints.

3.5. Inter-Temporal Costs

Inter-temporal costs have to be incorporated to ensure that transitions from one operating point to the next are feasible (i.e., managing flexibilities provided by the DER) and economical. These costs do apply independently to DER or assets while spanning multiple periods in particular, if multiple tariffs are defined along the horizon.

The costs for the controllable resources such as BESS and PVs follow a conditional operation regarding their mode of operation. The cost function, for instance, of charging a BESS is expressed by a linear convex function depending on the quantity of energy consumed—see the Figure 3a scenario when charging and discharging are equally priced. Any time step (e.g., charging or discharging) is dealt by imposing a Cost Constrained Variable (CCV) to represent the proper cost. The piecewise linear cost function $c_{BESS}\left(x\right)$—owned by the DSO—represented by the red line on Figure 3a—is substituted by an auxiliary variable $y_{BESS}$ and a set of linear constraints. These linear constraints form a convex region with the $c_{BESS}\left(x\right)$, setting the $y_{BESS}$ always to lie in the epigraph of the cost function. This auxiliary variable is onwards reflected to the objective function:

$$c_{BESS}\left(x\right)=\left\{\begin{array}{ccc}{\pi }_{BESS}·x& \mathbf{if}x\ge 0& \left(\mathrm{discharging}\right),\\ -{\pi }_{BESS}·x& \mathbf{else}x<0& \left(\mathrm{charging}\right),\end{array}\right.$$

(39)

where ${\pi }_{BESS}$ corresponds to price of utilizing the BESS either for charging or discharging. The necessary linear constraints for the auxiliary variable are the following:

$$\begin{array}{ccc}{y}_{BESS}& \ge & {\pi }_{BESS}·x,\\ y_{BESS}& \ge & -{\pi }_{BESS}·x.\end{array}$$

(40)

Accordingly, the cost functions and their subsequent auxiliary variables are defined. Note that, for EV or domestic BESS, their cost function can have an arbitrage regime, in the sense that charging has a profitable impact on the objective function while discharging is penalizing it—Figure 3d. Note that all Figure 3 are indicative for a particular period $\tau$; different functions might be regarded for different time steps implementing demand response schemes based on variable pricing schemes.

4. Case Study Synopsis

The network selected as a case study to perform the validation of the proposed scheme belongs to the IEEE benchmarked LV European network [51]. The network corresponds to a real British low voltage feeder connected to the MV grid through a transformer of nominal power 800 kVA (Figure 4). The transformer is modified and considered to be 20 kV to 400 V with nominal power of 150 kVA, since only one feeder is considered in the benchmark as well. The service cables to the 55 end-consumers are also included in the grid representation.

In all scenarios, a three-phase centralized BESS is assumed to be connected at node 101 (Figure 4). The BESS capacity is 90 kWh and the maximum charging and discharging rate is 45 kW. This BESS is assumed to have an initial SoC of 0.40 with $\underline{SoC}=0.1$ and, at the end of the optimization horizon, it has to be equal to the initial, $So{C}_{H\tau }$ = 0.40. The power factor of BESS${}_{101}$ is considered as unitary for all the simulated cases and its charging and discharging efficiency is ${\eta }_{ch},{\eta }_{dch}=$0.95, [52].

The load and the microgeneration profiles used correspond to daily data for a summer period, which are extrapolated from a realistic data pool provided for the benchmarked grid, which can be in [51]. All the consumers are single-phase and their phase connection is depicted in Figure 4. All the microgeneration units are considered as single-phase PV rooftop installations that are connected to the same phase as the respective residential user. The simulation day corresponds to a representative summer week day where typically peak loading conditions occur, where the load and generation profiles are illustrated in Figure 5a,b, accordingly.

Four EV models considered are based on four different EV models Nissan Leaf, Chevrolet Volt, BMW i3, Tesla S, which present different technical features regarding the Battery Capacity and charging power, in addition to their driving efficiency, which are considered from [53]. Their characteristics are listed in

Table 1. Electric Vehicles’ models and characteristics.

EV Model	Battery Capacity [kWh]	Charging Power [kW]	Driving Efficiency (km/kWh)	End-User Owner
Nissan Leaf	24	4	6.7	249.2/861.1/264.3/522.2
Chevrolet Volt	16	3.75	3.75	327.3/755.2/886.2/906.3/780.3/ 619.3/899.2/337.3/701.3
BMW i3	22	11	7.2	785.2/225.1/314.2/320.3/ 817.3/702.2/178.2/73.1/342.3
Tesla S	60	11	6.7	563.1/47.2/208.3/682.2/406.2/ 248.2/458.3/83.2/349.1/289.1

. All Tesla S and BMW i3 models are charged with an Efacec HomeCharger 7.4 kVA, while the rest of the EVs through an Efacec HomeCharger 3.4 kVA [54]. Therefore, the maximum charging power of each EV is limited and driven by the home charger used.

Concerning the EV usage, a routine is built to emulate credible scenarios, which are fed with public statistical data by [55]. This routine aims to capture patterns pertaining the EV usage upon different trip purposes as well as the trip duration (minutes) and length (km). The resulting data reflect a realistic response for the EV behavior during a day of operation, standing on the assertion that EVs charge exclusively at home. Nonetheless, the available statistical data do not explicitly provide information if the trip is from or to a destination. Therefore, the data as suggested in [56] are split into starting a trip and ending a trip, and it is assumed that a driver starts and ends every trip at home; hence, the cumulative data of trips for different days among different purposes are illustrated in Figure 6a. To assess the SoC change per each EV model used, averaged values for the purpose of each trip is used and correlated with each EV model’s driving efficiency (Table 1).

The following assumptions are also regarded for the EV:

• Initial SoC for all EV models is $So{C}_0$ = 0.5, which is meant to be the same at the end of the horizon $So{C}_{H\tau }=So{C}_0$.
• The charging efficiency and discharging—when V2G—efficiency are considered 85% for all EV models.

Two EV charging strategies are considered in this study:

• “Dumb” charging or uncontrolled charging where the EVs are not incorporated within the proposed operational scheme. Such uncontrolled charing profiles are extracted using the data in probability density function presented in Figure 6b.
• “Smart” charging, where the EV owner communicates relevant data (i.e., flexibility as defined above) regarding their commute and accordingly its availability to be charged according to the proposed tool. The V2G mode services enable the option to utilize the EV essentially for grid services. These constraints are automatically incorporated in the multiperiod-OPF scheme as the generalized set of Equations (13d)–(13e), whenever the availability of the EV allows it. The availability of the EV to charge is considered along the day during their idle periods (i.e., parked at the owner’s house).

Concerning the case where EVs follow the dumb charging, their charging occurs based on the distribution function given by [57]. The time departure and arrival as well as the daily distance traveled by each are randomly selected for a summer week day as presented in Figure 7.

According to the standard EN-50160, the 10 min mean r.m.s voltages shall not exceed the statutory limits during 95% of the week. Meanwhile, all 10 min mean r.m.s voltages shall not exceed the range of the Vn% + 10% and Vn% − 15% (which corresponds to 253 V and 195.5 V for most European grids). Given the fact that the proposed control scheme uses 30-min averaged data resolution, the voltage limits are set in [0.95, 1.05] p.u.values [58].

The percentage of PV and EV integration refers to the proportion of end users that own such units. For instance, 35% of EV penetration (i.e., 30 EVs), where the charging point of the EV are indicated in the last column of Table 1.

The use of DER is prioritized through the weighted terms $c_k\left(\tau \right)$, in the sense that the operational tool attempts to manage the flexibilities by addressing any voltage issues and respects the secondary transformer’s rated power with the controllable DER assigned with the more economical combination of $c_k\left(\tau \right)$. These $c_k\left(\tau \right)$ can be attributed with real operational cost values to reflect monetary values for the use of flexibility. For this study, the weighted terms for the use of flexibility for each type of DER derive a merit order scheme settled as $c_{BESS}<c_{EV}<c_{V2G}<c_{QR}<c_{APC}$, which present respectively the price of using the BESS, incorporating an EV in the coordinated charging, the use of V2G mode and finally the use of reactive and active power by the microgeneration. In this way, the tool prioritizes the use of the centralized BESS that is owned by the DSO, avoids excessive active power curtailments, and the presence of the EVs restrains the dispatch of reactive power by the PVs, which is rather not effective for addressing voltage issues in LV grids (i.e., $R>X$). Note that, in this study, the V2G operation is set slightly cheaper than any control of the microgeneration (i.e., APC or QR).

Prior to the presentation of the results within the proposed control approach, an exploration of discrete incremental integrations of PV and EV are given. For both cases, no control was deployed. In Figure 8a, the impact in voltage magnitudes is depicted; notable voltage drops can be observed for an EV integration about 35% (i.e., 20 EVs). Higher EV integration of 55–65% present severe voltage drops as well as the overloaded condition for the secondary substations (see Figure 8b). In Figure 9a, the analogous scenario for PV integration is presented. Overvoltages appear in scenarios with more than 35% of PV integration, where the reverse power flow towards the secondary substation can be also viewed through Figure 9b. One can notice that, in a mixed scenario with PVs and EVs, overvoltages will typically appear during morning and evening hours, whereas voltage sags will arise at late hours when most regular charging of EVs occur (see Figure 6a).

The explored scenarios in which the proposed control scheme is validated are defined in

Table 2. Scenarios’ description.

Scenario	Case 01 (C1)	Case 02 (C2)	Case 03 (C3)	Case 04 (C4)	Case 05 (C5)	Case 06 (C6)
EV [%]	0	0	0	35	55	65
PV [%]	55	73	85	0	0	55
BESS 101	✔	✔	✔	✔	✔	✔

. Analytical information regarding the points of connection for the micro-generation and EVs are given in Appendix A. The scenarios are selected in order to validate that the MACOPF could allow high integration of microgeneration avoiding overvoltages and overloading of the transformer; secondly, mixed scenarios of DER integration including EVs are assessed by the coordinated control amongst them.

5. Results

The assessment of the proposed control framework takes places in all the scenarios for a day-ahead deterministic planning of operation with a time-step of $\mathrm{\Delta }\tau =$ 30 min, (i.e., 48 time steps). The available active measures for operational purposes are active power control of the centralized three-phase BESS, the coordinated charging of the EVs in addition to V2G mode of operation where both are considered once the EV is available (i.e., parked at house premises), and the control active and reactive power of the microgeneration. The setup for the controller considers minimum power factor for microgeneration 0.9 and maximum allowed curtailment 15% of the injected power by each unit.

5.1. Cases C01–C03

In scenario C01, with 65% (36 PV units) overvoltages (up to 1.062 p.u) do arise along the grid due to notable active power injections. Additionally, reversed power flow is also present, increasing the loading conditions of the secondary transformer up to about 90%. In Figure 10a, the actions taken along the horizon period are illustrated. Among the decisions, BESS${}_{101}$ is essentially charged during sunny periods, reducing the reversed flows to the substation and restraining the voltages. Ultimately, the voltages issues are addressed in coordination with reactive power dispatch by the microgeneration (Figure 10a). Note that, for the same scenario, the use of the control approach solely with APC derives 5.3 kWh enabled, which corresponds to minimum curtailments since the controller is centralized (Figure 10b); hence, dealing with local P–V droop control would lead to higher curtailments since there is no topological confluence among PV and the decisions are merely based on the voltage at the point of connection of the inverter.

Scenario 02 refers to the point of PV integration where overvoltages and secondary substation are within admissible bounds without curtailing any active power by the PV units. The subsequent controls derived in this case refer to the coordinated operation of the centralized BESS with reactive power dispatch of a set of microgeneration units. The control setpoint of the BESS and its SoC for the day of operation are illustrated in Figure 11a–c.

Given the selected options for the controller (i.e., minimum PF = 0.9 and maximum APC 15%), the proposed control is capable of mitigating overvoltages and ensuring rated power for the secondary substation up to 85% of PV integration. Analytically, the actions taken for this scenario are illustrated in Figure 10d.

Table 3. Resulting curtailed active power and dispatched reactive power for the microgeneration units.

Scenario	Curtailed Energy [kWh]	Reactive Energy [kVArh]
Case 01	0	3.94
Case 02	0	7.43
Case 03	34.5	59.3

presents the anticipated curtailed power and dispatch reactive power within the day of operation for C01–C03.

Figure 12 presents all cases C01–C03 and the corresponding transformer loading conditions if the decision derived by the MACOPF are followed. It is noted that in all cases the MACOPF delivers actions which restrain the reversed power to the secondary substation. Only in scenario C03 is there a slight relaxation of the respective constraint since all of the active measures have been optimally utilized.

Excessive active power curtailment may indirectly mean high compensation costs for the DSO. The installation of a BESS can reduce the need of APC as shown in scenarios 01–03. In particular, scenario 01 presented that, solely based on APC, 5.3 kWh at least would be curtailed with local P–V droop based control. This fact may justify the investment from the DSO side when a high integration of PV prevails.

As a remark stands, the current analysis leans on averaged values of 30-min resolution; therefore, notable overvoltages may appear in higher resolution analysis for less PV integration. The latter justifies the tight limits posed for the voltage magnitudes for the purpose of this study. Concurrently, the importance of coordinated actions will be observable even in scenarios with less PV integration.

5.2. Cases C04–C06

In this part of the study, the EV integration is introduced within the LV grid. Initially, cases 04–05 examine the integration of 20 and 30 EV accordingly. The EV are randomly distributed along the end-consumers, only following the assumption that no more than one EV can be connected at each consumer. The point of connections for each scenario are given in Appendix A.

Case 04 presents voltage drops along the early hours (12:30 a.m. to 4:00 a.m.) but most significantly and below the posed lower limits in the late hours where the peak load is met (9:00 p.m.–10:00 p.m.). In this scenario, the MACOPF derives an optimal scheduling (Figure 13b) of the DSO’s BESS in order to address the undervoltages in the late hours. The coordinated charging for the current case is not deemed to be needed according to the proposed control. Nonetheless, it can be seen from Figure 13a that a slight increase of the peak load (i.e., hereby 6:30 p.m. to 11:00 p.m.) might trigger critical voltage sags and possible transformer loading conditions. It should be stressed that the 35% integration of EV, following dumb charging, leads to an increase of 22% of the peak load compared to the base case.

The following scenario C4 emulates a case with 55% of EV integration. In this case, the importance of multiple DER coordination is stressed since critical voltage sags can be noticed if EVs follow the dumb charging operation. The MACOPF yields a scheduling for the BESS${}_{101}$, which essentially supports during the late hours (i.e., 8:00 p.m.–9:30 p.m.) in addition to some EVs that are decided to operate in V2G mode. The peak load conditions along this interval are deteriorated due to the EVs that are added to consumption while returning from their trips. The coordinated charging during the early hours avoids the voltage drop along the distribution feeder, which would be noticed in the uncoordinated charging (Figure 14a).

The centralized BESS${}_{101}$ is not adequate to entirely mitigate voltage problems due to the fact that some of the EVs are connected to the ending point of the feeder. Therefore, V2G mode of operation (i.e., for EVs that are willing to provide it) and coordinated charging efficiently added to the coordinated operation to maintain the voltages within the bounds. The SoC for all EVs are illustrated in Figure 14b. The V2G mode of operation comprises generally an effective measure to address undervoltages in LV grids due to the fact that active power is purely injected ($R>X$). Nevertheless, proper spatial and temporal distribution of the available EVs to provide such service is necessary. The spatial distribution refers to the bus and the phase of connection of the EV, whereas the temporal distribution to the availability along the horizon.

The following scenario C06 presents a case with increased PV and EV integration (i.e., 55% and 65%, respectively). The presence of PV installation injects notable active power during the sunny period (10:00 a.m. to 3:00 p.m.), whereas, in the rest of the time slots, voltage drops are observed due to the dumb charging. Notice that, in the current study, the PVs and EVs are randomly distributed, in the sense that a PV installation is not necessarily connected at a point where an EV is placed (see

Table A1. SPoints of PV Connections for the Presented Scenarios.

Points of PV connections	522.2/388.1/178.2/ 676.2/639.2/337.3/ 701.3/614.3/562.1/ 682.2/70.1/556.3/ 629.1/47.2/349.1/ 563.1/264.3/458.3/ 249.2/289.1	611.1/74.1/320.3/ 73.1/276.2/225.1/ 327.3/387.1/619.3/ 702.2	702.2/502.1/342.3/ 208.3/539.3/ 688.2/406.2	248.2/83.2/314.2/ 896.1/785.2/900.1/ 899.2/755.2/780.3/ 898.1/813.2
Scenarios [% of PV]	30%
	55%
	65%
	85%

and

Table A2. Points of EV Connections for the Presented Scenarios.

Points of EV connections	327.3/835.3/785.2/563.1/ 755.2/249.2/225.1/47.2/ 886.2/898.1/314.2/208.3/ 906.1/861.1/320.3/682.2/ 780.3/406.2/817.1/248.2	619.3/860.1/702.2/ 458.3/899.2/264.3/ 178.2/83.2/337.3/ 556.3	73.1/349.1/701.3/ 522.2/342.3/289.1
Scenarios [% of EV]	30%
	55%
	65%

). In addition, the fact that the majority of EVs progress in the sunny hours essentially implies that only a few could be available to be shifted for charging in these time slots. The subsequent results for this case are illustrated in Figure 15a. In this scenario, as it can be seen from the SoC of EVs in Figure 15b, the EVs are charged within a coordinated way from 12:00 a.m. to 5:00 a.m. where light loading conditions prevail in the grid avoiding the undervoltages that are present within the dumb charging. Nonetheless, due to increased charging demand, idle EVs provide support at the slot 5:30 a.m. as well as at 5:00 p.m. to 6:00 p.m., maintaining the voltages within the bounds. Accordingly, during sunny periods, shifted EV charging can be noticed as well.

6. Conclusions

This work presents a tool that is capable of providing support to the DSO decision for the operation of LV three-phase distribution grids with increased integration of DER. The computational techniques proposed based on the explicit calculation of the first and second order derivatives (i.e., Jacobian and Hessian of the Lagrangian and the objective function), in addition to pivotal adjustments in the Jacobian, ensure a tractable optimal control based on the exact AC power flows.

The proposed centralized scheme ensures admissible voltage profiles by minimizing the active power curtailments of microgeneration through the efficient coordination of DER—maximizing, in this sense, the integration of microgeneration. The coordinated operation among the DER units reassures in the presented study up to 73% integration of microgeneration avoiding any curtailment brought by PV units. The EV integration can be also maximized if coordinated charging is adopted within the MACOPF, which essentially can ensure admissible voltages and normal loading condition for the transformer. The V2G mode of operation can be regarded as important when a high integration of EV takes place with increased peak load conditions in the grid.

For all simulated cases, one can conclude that the radial configuration of LV networks present overvoltage issues particularly in buses—electrically—furthest from the substation since the high resistance of the lines leads to the aggravation of them. In addition, distant nodes at the ending point might face significant voltage drops (i.e., if EVs are present) or even overvoltages (i.e., if PVs are present). Therefore, the placement of a BESS in proper location (i.e., adjacent to nodes with voltage issues) along the grid is rather crucial since the installation adjacent to the secondary substation provides mainly reduced loading for the transformer rather than voltage support to the furthest points.

The proposed control scheme is based on deterministic analysis for the planning in a day-ahead scale for the operation of the grid. Meanwhile, the scheme can be deployed only by the subsequent communication technologies (for online implementation), together with forecasted data and power flow-state estimation tools (i.e., if topology is not known).