Optimal Allocation of ESSs in Active Distribution Networks to achieve their Dispatchability
Ji Hyun Yi Rachid Cherkaoui Mario Paolone

	Non-standard interactions versus planet-scale neutrino oscillations
    Yi-Xing Zhou
    December 30, 2023
===================================================================




This paper presents a method for the optimal siting and sizing of energy storage systems (ESSs) in active distribution networks (ADNs) to achieve their dispatchability. The problem formulation accounts for the uncertainty inherent to the stochastic nature of distributed energy sources and loads. Thanks to the operation of ESSs, the main optimization objective is to minimize the dispatch error, which accounts for the mismatch between the realization and prediction of the power profile at the ADN connecting point to the upper layer grid, while respecting the grid voltages and ampacity constraints. The proposed formulation relies on the so-called Augmented Relaxed Optimal Power Flow (AR-OPF) method: it expresses a convex full AC optimal power flow, which is proven to provide a global optimal and exact solution in the case of radial power grids. The AR-OPF is coupled with the proposed dispatching control resulting in a two-level optimization problem. In the first block, the site and size of the ESSs are decided along with the level of dispatchability that the ADN can achieve. Then, in the second block, the adequacy of the ESS allocations and the feasibility of the grid operating points are verified over operating scenarios using the Benders decomposition technique. Consequently, the optimal size and site of the ESSs are adjusted. To validate the proposed method, simulations are conducted on a real Swiss ADN hosting a large amount of stochastic Photovoltaic (PV) generation.



Energy storage systems, optimal power flow, active distribution networks, resource planning, dispatchability











§ INTRODUCTION
















Power balancing is becoming an increasingly challenging task due to growing volatility in power systems introduced by the high penetration penetration of non-dispatchable and stochastic generation.
 The spinning reserve requirement is growing to tackle not only the peak demand but also the stochasticity of prosumption, which is defined as the load consumption minus the generated power, causing significant cost increases for both transmission system operators (TSOs) and distribution system operators (DSOs). In the distribution grid, the prosumption forecast uncertainties are originated not only by the load consumption but also by the distributed RES. Thus, the deviation of the network infeed from the prescheduled power intake from the grid, called dispatch error, should be compensated by the system operating reserve, exposing DSOs to high imbalance cost. 
There has been increasing interest in using ESSs as a flexible resource to compensate for the dispatch error, thereby enhancing the dispatchability of the ADNs and mitigate the aforementioned financial risk. The dispatchability identifies the capability of a resource, or a network, to control the realized active power flow through the resource, or the network, to follow a pre-defined power schedule with high fidelity. Thanks to the ADN dispatchability, the responsible DSO ensures its profit by mitigating imbalance penalties <cit.>. In this respect, the optimal control strategies of ESSs to improve the dispatchability of the stochastic resource <cit.>, or distribution feeder <cit.> has been addressed. 

However, the efficient deployment of ESSs significantly depends on the investment planning. The work in <cit.> proposed sizing of a wind farm ESS to achieve its dispatchability. The ESS dispatch strategy is coupled with the assessment of its capacity and expected lifetime based on the confidence level of the power output w.r.t the predefined schedule. To ensure the economic optimality of decision-making, it is worth investigating the economic value of ESSs investments by taking into account their operational benefits and installation costs. The work in <cit.> proposed a method for the optimal allocation of ESSs through a cost-benefit analysis considering mitigating the distribution utilities' risk for compensating the gap between the actual prosumption and the purchased energy from the forward market. <cit.> is an extension of <cit.>, where different agents, such as DSOs, wind farms, solar power stations and demand aggregators in ADNs, determine the capacity of ESSs based on the game theory approach to reduce the transaction costs risk due to the resources' forecast errors.

It is also worth observing that the capacity and placement of ESSs should comply with the characteristics of the chosen control strategy and the subsequent operational conditions of the system. In this regard, the optimal power flow (OPF) that accurately models operation and control of the ADN, should be embedded in the ESS planning tool. However, OPF-based ESSs planning is inherently burdensome to solve due to its non-convexity. There exists meta-heuristic methods such as fuzzy PSO algorithm used in <cit.> and <cit.>. However, these solution approaches cannot guarantee the global optimum or even a feasible solution. To tackle the non-convexity of the OPF problem, convexification approaches, such as the semi-definite programming (SDP) <cit.> or the second-order cone programming (SOCP) relaxation <cit.>, have emerged as a solid and rigorous alternative. 
The SOCP relaxation proposed in <cit.> has been implemented for the optimal ESS allocation in radial grids due to its superior computational efficiency referring to the SDP relaxation. In <cit.>, an ESS planning strategy was developed relying on the SOCP-OPF model with the objective of providing ancillary services to the TSO, and to cope with the wind variability. Authors of <cit.> further utilized the SOCP relaxed model for the ESSs allocation and operation to minimize the grid losses and imported power in the ADN.  
The work in <cit.> tackled the ESS planning and operation problem by decomposing it into two stages: first stage determines the total ESS size to prevent grid constraints violations due to PV power imbalance, and the second stage allocates ESSs with optimal sizes by employing SOCP-OPF to minimize the energy cost.  

Meanwhile, the Authors addressed that the objective function and constraints should satisfy some necessary conditions to guarantee the exactness of the SOCP-OPF solution, and suggested a formula to verify the exactness a-posteriori. The main drawbacks were explicitly underlined in <cit.> by the fact that the exactness of the solution cannot be guaranteed especially in the presence of reverse line power flows and for the cases where the upper bound of nodal voltage and the line ampacity constraints are binding.This brings significant limitations on the applicability of the method to ADNs hosting distributed generation units with large capacities. Moreover, the model neglects the transverse elements of the lines, which can bring an infeasibility of the solution especially when ADNs are composed of coaxial underground cables. 
The work in <cit.> solves this problem and proposed the Augmented Relaxed OPF (AR-OPF) to convexify the AC-OPF for radial grids. Their contribution demonstrates that, in the AR-OPF problem comprising an objective function strictly increasing with the grid losses, the conditions for the exactness of the solution are mild and hold for realistic distribution networks. The AR-OPF was implemented in the subsequent works on the optimal ESS planning problem while embedding grid reconfiguration with the objective of minimizing grid losses, voltage deviation and the line congestion <cit.>.  

In this paper, we propose an operation-driven planning strategy of ESS to achieve the dispatchability of the ADN based on the AR-OPF model. The objective of achieving disatchability requires a substantial and non-trivial modification on the AR-OPF as well as on the solution approach proposed in <cit.> in order to reach the exactness of the relaxed OPF. Moreover, we formulated the sizing problem into two blocks by modifying the objective term, constraints and variables related to the dispatch error. Meanwhile, we apply the Bender’s decomposition to handle the multi-layered decisions with numerous scenarios <cit.>. The contributions of the paper are two as follows.

	

  * The optimal allocation of ESS is determined based on an exact convex model of the OPF to address the dispatchability of the ADN in the presence of prosumption uncertainty, while accurately reflecting the operational condition of the ADN.

  * The structure of the planning problem and the objective functions are formulated accounting for the necessary conditions to guarantee the optimality and the exactness of the OPF relaxation.


The paper is organised as follows: in Section II, we introduce the structure of the optimization problem and explain the key parts in detail. In Section III, the proposed problem formulation and solution approach are described. Section IV contains a detailed application example referring to the planning of ESSs into a real ADN. Finally, Section V concludes the paper.


§ PROBLEM STRUCTURE





 §.§ System description

Buses excluding the slack bus in an ADN are indexed with l ∈ℒ. Lines whose downstream bus is bus l is are also indexed with l ∈ℒ. In this paper, optimal allocation of ESS is determined based on the operation of the ADN over the planning horizon Y, assuming that there is no change in system conditions over years. Each year, we classify days into day-types indexed with d ∈𝒟. The uncertainty in prosumption in each day-type is represented by operating scenarios indexed with ϕ∈Φ_d with probability λ_ϕ d for each scenario. The dispatch interval is identified by the index t∈{1,...,T}=𝒯, where T corresponds to the scheduling horizon of the problem. Time indices are separated by a constant timestep Δ t. The dispatch problem at day d accounts only for the active power. Each node of the ADN has a non-dispatchable aggregated prosumption (s_lϕ t) defined at each scenario and time interval. The aggregated ADN active powers through the grid connecting point (GCP) to an upper layer grid in all scenarios (P_1ϕ t, ∀ϕ∈Φ_d) are expected to follow a day-ahead determined daily dispatch plan at each time interval (DP_td) derived with the support of a forecasting tool. Then, at each node where an ESS is allocated (i.e., U_l=1, where U_l∈{0, 1} is installation status of the ESS at bus l), it is dispatched at each scenario and time interval according to active power (p^E_lϕ t) and reactive power (q^E_l ϕ t). The dispatched active power compensates for the gap between active power infeed and the dispatch plan, which results from the deviation of realized prosumption from the prosumption prediction. Consequently, the observed dispatch error at GCP can be minimized.

 
In summary, the ESS allocation problem to achieve the dispatchability of the ADN under study is a two-stage decision process: the first stage which deals with the binary decision variables on the location of the ESS (U_l) and the continuous decision variables on the capacity of the ESSs energy reservoirs (C_l) and their power rating (R_l), and the second stage which deals with daily dispatch problems, determining the decision variables on the ESSs active and reactive power for all operating scenarios. 


 §.§ Dispatch plan of the distribution feeder


The operational benefit of ESSs allocation for the ADN dispatchability is evaluated with sets of operationg scenarios, where each set refers to a typical day. The set of day-types is supposed to be pre-selected by the modeler.  

The advantage of using scenarios is related to the use of generic parametric and non-parametric distributions of prosumptions. In this paper, the active and reactive prosumption scenarios for each day are generated with the assumption that the prosumption profile follow a normal distribution. Therefore, the average of the prosumption over the scenarios is equal to the given prosumption prediction, and the losses predictions are calculated by averaging the losses over the scenarios <cit.>. It is understood that the variables with subscript l, ϕ, t are defined for l∈ℒ, ϕ∈Φ_d, t∈𝒯, respectively. p_lϕ t is aggregated active prosumption at bus l, scenario ϕ, time t, and f_lϕ t is square of current magnitude causing losses in line l, scenario ϕ, and time t. r_l is the total longitudinal resistance of line l. p̃_ltd and f̃_ltd are avereage of the active load and squared current causing losses over all scenarios at bus l, time t, and day d, respectively. Δ p_lϕ t and Δ f_lϕ t are deviation (we would call it as error hereafter) of prosumption and squared current causing losses for scenario ϕ and time t from p̃_ltd and f̃_ltd, respectively. eq:prosumption prediction and eq:losses prediction express the prosumption at each bus l and the line losses of the line l by these variables.   
In eq:dispatch plan, a daily dispatch plan DP_td follows the predicted point of the total prosumption considering the losses.

    p_l ϕ t=p̃_ltd-Δ p_l ϕ t, ∀ l,∀ϕ, ∀ t 
        r_lf_l ϕ t=r_lf̃_ltd-r_lΔ f_l ϕ t, ∀ l,∀ϕ, ∀ t  
        DP_td=∑_l∈ℒ(p̃_ltd+r_lf̃_ltd), ∀ t

In this context, the dispatch error with no dispatchable resources in ADN is formally defined as the total error of prosumption plus the line losses over the buses/lines as shown in the left-hand side of eq:error expression. The installed ESSs in ADN can be dispatched with acive power (p_l ϕ t^E) to compensate for the error. 

 ϵ_l ϕ t represents the residual dispatch error that cannot be covered at bus l, scenario ϕ, and time t. To quantify the covered or not covered error by the ESSs, we say that the sum of the two parts, shown in the right-hand side of eq:error expression, is equal to the dispatch error in case of no dispatchable resources.  

    ∑_l∈ℒ(Δp_l ϕ t+r_lΔ f_l ϕ t)=∑_l∈ℒ(ϵ_l ϕ t+p_l ϕ t^E), ∀ϕ, ∀ t

Finally, the objective term to minimize the power deviation from the daily dispatch plan for all scenarios and time intervals is expressed as follows.

    _p^E∑_t∈𝒯∑_ϕ∈Φ_dλ_ϕ d|∑_l∈ℒϵ_l ϕ t|

Once the OPF problem including the constraints related to the dispatchability is solved, we can calculate the capability of ADN with the allocated ESS sizes to reduce the dispatch error. In this context, we introduce a dispatchability index called leftover dispatch error rate (LDER), and defined as θ_ϕ td. It represents the ratio between the resulting dispatch error and the anticipated dispatch error in case of no ESS at scenario ϕ and time t for the daily operation on day d, as expressed in eq:error reduction rate_def (* here indicates that it is the identified optimal solution). In the denominator, indicating the dispatch error without ESS, the error regarding the grid losses is ignored because its magnitude is negligible compared to that of the prosumption error. 

    θ_ϕ td=|∑_l∈ℒϵ_l ϕ t^*|/|∑_l∈ℒΔ p_l ϕ t|,∀ϕ, ∀ t
  


 §.§ Augmented Relaxed Optimal power flow


For a radial power network, the power flow equations are given in eq:complex power balance up-eq:relaxed internal current. The variables with subscript l are defined for l ∈ℒ. The upstream bus of bus l is notated as up(l). 𝐆 is the adjacency matrix of the network, where 𝐆_k,l is defined for k,l∈ℒ and 𝐆_k,l=1 if k=up(l) or 0 if not. S_l^t=P_l^t+jQ_l^t indicates the complex power injected into line l from node up(l). S_l^b=P_l^b+jQ_b^t is the complex power injected into node l from line l. s_l=p_l+jq_l represents the given prosumption at node l. v_up(l) and v_l represents the nodal voltages at node up(l) and l, respectively, whereas f_l is the square of the longitudinal (i.e., internal) current through line l which produces grid losses through the line longitudinal impedance z_l=r_l+jx_l. b_l indicates the shunt susceptance of line l. S_l^t and S_l^b are determined by eq:complex power balance up, and eq:complex power balance down, respectively. eq:nodal voltage determines the magnitude of nodal voltages, where ℜ(.) represents the real parts of a complex number. The current that produces losses at line l is originally defined to be equal to the right-hand side of eq:relaxed internal current. However, the equality is replaced with an inequality, as shown in eq:relaxed internal current by applying the SOCP relaxation <cit.>. It is noteworthy that the equality of eq:relaxed internal current only holds when the terms in the objective function are strictly increasing with respect to the internal current f_l.

    S_l^t=s_l+∑_m∈ℒ𝐆_l,mS_l^t+z_l f_l-j(v_up(l)+v_l)b_l, ∀ł
        S_l^b=s_l+∑_m∈ℒ𝐆_l,mS_l^t, ∀ł
    
    
        v_l=v_up(l)-2ℜ(z_l^*(S_l^t+j v_up(l)b_l ))+|z_l|^2f_l, ∀ł
    
    
        f_lv_up(l) ≥ |S_l^t+j v_up(l)b_l|^2, ∀ł


In order to avoid any inexact solution (i.e., any solution that makes the left-hand side of eq:relaxed internal current strictly greater than the right-hand side, and thus without physical meaning), the Authors of <cit.> introduced auxiliary variables to formulate the AR-OPF model. These are f_l, S_l=P_l+jQ_l, v_l and Ŝ_l=P̂_l+jQ̂_l, which stand for the upper bounds for internal current, complex power, voltage and lower bounds for the complex power, respectively. v^max and v^min are upper bounds and lower bounds of the squared nodal voltage magnitude, respectively. I^max_l indicates upper limit on the squared current of line l. P^max_l and Q^max_l are upper limits of active and reacitve power flows for line l, respectively.
 The branch power flow, nodal voltage, and current equations are defined as well with the set of the auxiliary variables, as shown in eq:lowerbound complex power balance up-eq:upperbound internal current2. eq:lowerbound complex power balance up and eq:lowerbound complex power balance down indicate the lower bound of branch power flow at the sending end and the receiving end of line l. The upper bound nodal voltage is determined correspondingly with eq:nodal upperbound voltage. Likewise, the branch power flow equations for upper bound power flow variables are shown as eq:upperbound complex power balance up and eq:upperbound complex power balance down. eq:upperbound internal current and eq:upperbound internal current2 express that the upper bound of the current f_l should be decided by the maximum of absolute complex power flow from both sides of line l. The voltage constraint, the ampacity constraint from the sending end and the receiving end are modeled as in eq:voltage constraint-eq:ampacity from downside. Finally, eq:upperbound power range are added to complete the set of equations required to guarantee the exactness. 

    Ŝ_l^t=s_l+∑_m∈ℒ𝐆_l,mŜ_l^t-j(v_up(l)+v_l)b_l, ∀ł
    
    
    
    
    
    
    
    
       Ŝ_l^b=s_l+∑_m∈ℒ𝐆_l,mŜ_l^t, ∀ł
    
    
       v_l=v_up(l)-2ℜ(z_l^*(Ŝ_l^t+jv_up(l)b_l)), ∀ł
    
    
       S_l^t=s_l+∑_m∈ℒ𝐆_l,mS_l^t+z_lf_l-j(v_up(l)+v_l)b_l, ∀ł
    
    
    
    
    
    
    
    
       S_l^b=s_l+∑_m∈ℒ𝐆_l,mS_l^t, ∀ł
    
    
       f_lv_l    ≥|max{|P̂_l^b|,|P_l^b|}|^2
       +|max{|Q̂_l^b-j v_l b_l|,|Q_l^b-j v_lb_l|}|^2, ∀ł
    
    
       f_lv_up(l)   ≥|max{|P̂_l^t|,|P_l^t|}|^2
       +|max{|Q̂_l^t+jv_up(l) b_l|,|Q_l^t+j v_up(l)b_l|}|^2, ∀ł
    
    
    
    
        v^min ≤ v_l,  v_l ≤ v^max, ∀ł
    
    
        I_l^maxv_up(l)≥|max{|P̂_l^t|,|P_l^t|}|^2+|max{|Q̂_l^t|,|Q_l^t|}|^2, ∀ł
    
    
    
    
    
    
    
    
    
    
    
        I_l^maxv_l≥|max{|P̂_l^b|,|P_l^b|}|^2+|max{|Q̂_l^b|,|Q_l^b|}|^2, ∀ł
    
    
       P_l^t ≤ P_l^max,     Q_l^t ≤ Q_l^max, ∀ł

For the sake of readability, the equations mentioned above are grouped and represented by Θ(φ,κ)≥0 where φ={S^t,v,f,Ŝ^t,v,f, S^t,s} is the set of variables and κ={𝐆,z,b,v^max,v^min,I^max, P^max, Q^max} is the set of parameters. The notation without subscript corresponds to the vector of variables and parameters for all buses/lines. 

As already quantified in <cit.>, the set of grid constraints employing the auxiliary variables slightly shrinks the original feasible solution space, removing solution space related to undesirable or extreme operation points of the network near the upper bound of nodal voltage or lines' ampacity limits.[We show that the compression of the solution space caused by the augmented constraints is small by following the steps of numerical analysis reported in <cit.>. Under the case study shown in sec:simulation, We make one of two operating constraints (voltage upper bound constraint and ampacity constraint) binding at one node or line and relax the other one to find the difference between the physical state variables (nodal voltage-magnitudes and original current flow) and corresponding auxilary variables. The difference between the nodal voltage-magnitude and the auxillary one is 0.001%. The difference between the original current flow and the auxillary one is equal to 0.2%.]













 §.§ Energy storage systems


The ESS investment is modeled through eq:available p rating-eq:c-rate.  
 𝒞^max_l/𝒞^min_l is maximum/minimum possible ESS energy reservoir capacity at bus l. ℛ^max_l and ℛ^min_l represent maximum and minimum possible ESS power rating capacity at bus l. In reality, available power ratings and energy capacities are often restrained as seen in eq:available p rating and eq:available e capacity due to various physical constraints involving, for instance, manufactural or geographical factors. CR_min is the minimum value for the rate at which ESS is discharged relatively to its maximum energy capacity. The power rating and energy reservoir is determined considering this relationship as shown in eq:c-rate. 

    ℛ^min_lU_l≤ R_l≤R^max_lU_l, ∀ł
    
    
       𝒞_l^minU_l≤ C_l≤𝒞_l^maxU_l, ∀ł
    
    
        R_l ≤C_l/CR_min, ∀ł

The operational characteristic of ideal ESS is modeled with equations eq:power ESS-eq:SOE limit. The variables with subscript l are defined for l ∈ℒ. E_lt^E is state-of-energy (SoE) of ESS installed at bus l for time t. E^max_l/E^min_l is maximum/minimum allowed SoE level. The capability curve of a given ESS is linearized by constructing an inscribed square within the original circular capability curve defined by the maximum complex power of the ESS (see eq:power ESS). The SoE level of ESS changes with the charge/discharge of the ESS at every time interval, as described in eq:dynamic SOE. Also, eq:SOE limit indicates that the SoE of the ESS should be within the SoE limits. 

    -R_l/√(2)≤ p_lt^E≤R_l/√(2),   -R_l/√(2)≤ q_lt^E≤R_l/√(2), ∀ l, ∀ t
    
    
    
        E_l(t+1)^E=E_lt^E+Δ t*p_lt^E, ∀ l, ∀ t
    
    
    
        E^min_l C_l ≤ E_lt^E≤ E^max_l C_l, ∀ l, ∀ t

In the interest of brevity, (<ref>)-(<ref>) are indicated by Ξ(η_t,ξ)≥0, ∀ t ∈𝒯 where η={p^E,q^E,E^E, U, R, C} is the set of variables and ξ={Δ t, E^min, E^max} is the set of parameters. The notation without subscript corresponds to the vectors of variables and parameters for all buses/lines.


§ PROBLEM FORMULATION

The objective of the problem is to determine the optimal sizes and sites of ESSs so that the active power through the GCP follows the dispatch plan with minimal deviation. However, it is clear that the dispatch error described as eq:modified dispatchability in sec:achieving dispatchability does not increase while the total grid losses increase (a necessary condition of the AR-OPF). Therefore, the exactness of the solution cannot be guaranteed if the objective value that corresponds to the objective term eq:modified dispatchability is significant in magnitude compared to objective term regarding the total grid losses in the objective function of the AR-OPF model.


Therefore, we propose to decompose the problem into two blocks each consisting of an OPF problem. In this way, we can exclude eq:modified dispatchability from the AR-OPF problem and convey it to another, approximated, OPF problem (the so-called 1st block problem), which aims to find the optimal level of dispatchability based on the ESSs investment cost and the imbalance penalty.

 The determined dispatchability level, defined as LDER, is then introduced to an AR-OPF model-based problem (the so-called 2nd block problem) as an index with which the ADN has to comply with.
 

The whole algorithm of the proposed approach is illustrated in fig:full_algorithm. In the 1st block problem, the optimal ESSs allocation, the daily dispatch plans, and the corresponding LDER are calculated employing linear approximated OPF ignoring the grid losses. Only the nodal voltage constraints are considered regarding the operational constraints, ignoring the ampacity limits to reduce the computational burden. Afterward, the outputs of the 1st block, which are the ESS allocation and the LDER, are used as inputs for the 2nd block problem.  

In the 2nd block problem, the objective is to refine the optimal allocation of the ESSs, considering several operating scenarios and achieving the same level of LDER calculated in the 1st block problem. In this respect, LDER is implemented as an additional constraint to an AR-OPF model, which considers the full AC-OPF as well as voltage constraints and line ampacity limits. Then, the size and site of the ESSs are iteratively adjusted thanks to the Benders decomposition technique. This iterative process starts initially with a feasibility check of the ESS allocation resulting from the 1st block problem. 



 §.§ 1st block problem

We minimize the investment and total penalty costs over the planning horizon to find out the optimal allocation of the ESSs and the optimal dispatchability level. We embed an approximated OPF constraints for all operating scenarios into a two-stage mixed-integer linear programming problem (see sec:system). The OPF is formulated by the linear Distflow model in which shunt elements are considered, whereas the grid losses are neglected. In this way, the reactive power generated by the shunt impedance of the lines is accounted for in the nodal voltage constraints. Meanwhile, neglecting the formulation of the squared current (f_l) (i.e., losses) is less likely to affect the feasible solutions in this stage since the ampacity constraint is ignored.

The dispatch plan follows the prosumption prediction, as shown in eq:distflow dispatch plan, while the prosumption deviation at node l is expressed as shown in eq:dist flow error expression. We substitute eq:prosumption prediction into the active power balance equation, resulting in eq:distflow power balance up dispatch. The lossless Distflow power flow at both sides of line l including the ESS power are expressed via eq:distflow power balance up dispatch-eq:distflow reactive power balance down. eq:distflow nodal voltage calculates the nodal voltage, which is governed by voltage constraint as shown in eq:distflow voltage constraint. All the variables within eq:dist flow error expression-eq:distflow voltage constraint are defined for ł∈ℒ, ϕ∈Φ_d, t∈𝒯, and d∈𝒟.


    DP_td=∑_l∈ℒp̃_ltd, ∀ t, ∀ d
    
       ∑_l∈ LΔ p_l=∑_l∈ℒ(ϵ_l+p_l^E), ∀ϕ, ∀ t, ∀ d 
    
        P_l^t=P_l^b=p̃_l-Δ p_l+p_l^E+∑_m∈ℒ𝐆_l,m P_l^t,∀ l,∀ϕ, ∀ t, ∀ d
    
        Q_l^t=q_l+q_l^E+∑_m∈ℒ𝐆_l,m Q_l^t-(v_up(l)+v_l)b_l,∀ l,∀ϕ, ∀ t, ∀ d
    
        Q_l^b=q_l+q_l^E+∑_m∈ℒ𝐆_l,mQ_l^t,∀ l,∀ϕ, ∀ t, ∀ d 
    
        v_l=v_up(l)-2ℜ(z_l^*(S_l^t+j v_up(l)b_l)),∀ l,∀ϕ, ∀ t, ∀ d
    
        v^min ≤ v_l ≤ v^max,∀ l,∀ϕ, ∀ t, ∀ d

 The objective function is defined as to minimize the investment cost (the first line of eq:objective function)  of ESSs and the annual penalty cost (the second line of eq:objective function) regarding the uncovered dispatch error over the planning horizon Y. The penalty cost of day d is the uncovered dispatch error over the operating scenarios for day d multiplied by ω_d, which is the cost coefficient for the imbalance. N_d stands for the number of days in a year that have the same prosumption profile as that of typical day d. Ω_1 and Ω_2 represents the set of control variables in the first and second stage decision process, respectively. The constraints regarding the ESS allocation and operation explained in sec:ESS are included (i.e., eq:available p rating-eq:c-rate, eq:ESS_constraints_1st), along with the linear approximated lossless OPF constraints (eq:distflow dispatch plan-eq:distflow voltage constraint).
By solving the problem, the optimal allocation of ESSs along with the dispatch plan is obtained, and the daily LDERs are calculated by eq:error reduction rate_def concerning all days with the index d ∈𝒟.

    _∀ U,C,R ∈Ω_1;
    ∀ S^t,v,s^E∈Ω_2   ∑_l∈ℒ (ℐ_c U_l+ℐ_p R_l+ℐ_e C_l)
       +Y∑_d∈𝒟N_dw_d∑_t∈𝒯∑_ϕ∈Φ_dλ_ϕ d|∑_l∈ℒϵ_l ϕ t|
 
subject to: eq:available p rating-eq:c-rate, eq:distflow dispatch plan-eq:distflow voltage constraint,  

    Ξ(η_ϕ t d,ξ)≥0,  ∀ϕ, ∀ t, ∀ d





 §.§ 2nd block problem

The dispatchability level is incorporated into the 2nd block problem as a constraint governed by the dispatchability index LDER. The objective of the 2nd block problem is to adjust the ESS allocation from the 1st block problem to the optimal site and size that can minimize the grid losses and unserved load. The system condition during the operation horizon is evaluated through solving the AR-OPF problem. Therefore, the 2nd block problem is formulated as a Mixed-integer second-order cone programming (MISOCP) problem. A convex SOCP model can be obtained from the non-convex MISOCP problem by using solution approaches that tackle the integer variables, such as the Branch and bound algorithm. However, the solution approach becomes computationally burdensome with the increase of integer variables. In this respect, we apply the Benders decomposition technique to decompose the 2nd block problem into a master problem and several parallel subproblems that each represents a daily OPF problem. The master problem determines the ESSs allocation, followed by the fitness evaluation of the determined allocations in the subproblems in terms of grid losses and unserved load. The unserved load takes values to ensure the feasibility of the subproblem regardless of the ESS allocation.

The initial step of the 2nd block is to solve parallel subproblems, checking the operational condition of ADN under the ESS allocation and the corresponding LDER calculated from the 1st block. The objective of the subproblem and the dual values of ESS allocation are computed and sent to the master problem to construct the first Bender's cut. Through multiple iterations between the master problems and the subproblems, the convergence is reached when the gap between the lower bound of the total cost (LB) and the upper bound of the total cost (UB) becomes less than a tolerance (i.e., 0.01% of UB). LB is determined from the master problem, whereas UB is determined after solving subproblems.



  §.§.§ Master problem

The formulation of the master problem is given in eq:objective function_of_master. The master problem computes the lower bound of the planning problem by summing the investment cost and the lower approximation of the subsequent expected subproblem costs. Each α_d represents the subproblem cost for days classified into each day type. It is initially bounded by α, which is the parameter given as the lower bound for the subproblem cost. n∈𝒩 is the index of benders iterations. In every nth iteration, Benders cuts represented by Γ_d^(n) (see eq:bender's_cut) are added as additional constraints for days d∈𝒟, as shown in eq:optimality_cut. The lowerbound of the total cost, so-called LB, is determined by the optimal solution of the master problem (i.e., LB=ℐ_c U_l^*+ℐ_p R_l^*+ℐ_e C_l^*+∑_d∈𝒟α_d^*).

    _∀, U,C,R,αℐ_c U_l+ℐ_p R_l+ℐ_e C_l+∑_d∈𝒟α_d


subjected to : (<ref>)-eq:c-rate,



    α_d ≥α, α_d ≥Γ_d^(n),  ∀ d, ∀ n



  §.§.§ Subproblem

In the subproblem associated with day d, a daily AR-OPF model with the time-step discretization of 15 minutes evaluates the operational advantages of ESSs while considering real operational conditions. The variables with subscript l, ϕ, t are defined for l ∈ℒ, ϕ∈Φ_d, t ∈𝒯, respectively.
The sufficiency of the ESS allocation is assessed by checking if the uncovered dispatch error (see eq:error expression) satisfies the LDER for day d as shown in eq:error reduction rate. As indicated in eq:active prosumption with unserved load and eq:reactive prosumption with unserved load, we introduce positive and negative unserved active load terms (ulp_l ϕ t^+, ulp_l ϕ t^-∈ℝ^+) and positive and negative unserved reactive load terms (ulq_l ϕ t^+, ulq_l ϕ t^-∈ℝ^+) to the active prosumption and reactive prosumption at bus l, scenario ϕ and time t, respectively. They correspond to the amount that should be curtailed from the prosumption to primarily comply with the LDER constraint along with other operational constraints, even in the case of insufficient capacity of ESSs. 

    |∑_l∈ℒϵ_l ϕ t|≤θ_ϕ t d|∑_l∈ℒΔ p_l ϕ t|,  ∀ϕ, ∀ t
        
        p'_l ϕ t=p_l ϕ t+ulp_l ϕ t^+-ulp_l ϕ t^-,  ∀ł, ∀ϕ, ∀ t
    
        q'_l ϕ t=q_l ϕ t+ulq_l ϕ t^+-ulq_l ϕ t^-, ∀ł, ∀ϕ, ∀ t
    
        p'_l ϕ t=p̃'̃_l ϕ t-Δ p'_l ϕ t,  ∀ł,∀ϕ, ∀ t
    
        DP_td=∑_l∈ℒ(p̃_̃l̃t̃d̃'̃+r_lf̃_ltd),  ∀ t 
    
       ∑_l∈ℒ(Δp_l ϕ t'+r_lΔ f_l ϕ t)=∑_l∈ℒ(ϵ_l ϕ t+p_l ϕ t^E),  ∀ϕ, ∀ t

The AR-OPF problem embedding the dispatchability for the subproblem is formulated by replacing eq:prosumption prediction with eq:prosumption prediction2. In other words, we replace the prosumption p_l by p'_l in the relevant equations. p̃'̃_̃l̃ is employed in place of p̃_̃l̃ to determine the daily dispatch plan, as shown in eq:dispatch plan2. Δ p_l of eq:error expression is substituted with Δ p'_l to build eq:error expression2. 
Similarly, p'_l and q'_l replace p_l and q_l in the active power balance equations formulated with the state variables (i.e., eq:complex power balance up, eq:complex power balance down) and the auxiliary variables (i.e., eq:lowerbound complex power balance up, eq:lowerbound complex power balance down, eq:upperbound complex power balance up, eq:upperbound complex power balance down).
The equations of AR-OPF model with ESSs are re-defined as eq:AROPF_ESS_constraints, where φ'={S^t,v,f,Ŝ^t,v,f, S^t,s'=(p'+p^E)+j(q'+q^E), ulp^+, ulp^-, ulq^+,
 ulq^-} is the set of variables and κ={𝐆,z,b,v^max,v^min,
I^max, P^max, Q^max} is the set of parameters. The ESS power is governed by the set of operational constraints as eq:ESS_operation.

    Θ'(φ_ϕ t',κ)≥0, ∀ϕ, ∀ t 
       Ξ(η_ϕ t,ξ)≥0,  ∀ϕ, ∀ t.
  


    P_l^t=p̃'̃_l+r_lf̃_l-Δ p'_l-r_lΔ f_l+p_l^E+∑_m∈ℒ𝐆_l,m P_l^t


    P_l^b=p̃'̃_l-Δ p'_l+p_l^E+∑_m∈ℒ𝐆_l,mP_l^t


    P̂_l^t=p̃'̃_l-Δ p'_l+∑_m∈ℒ𝐆_l,m P̂_l^t

 
    P̂_l^b=p̃'̃_l-Δ p'_l+p_l^E+∑_m∈ℒ𝐆_l,mP̂_l^t

 
    P_l^t=p̃'̃_l+r_lf_l-Δ p'_l-r_lΔf_l+p_l^E+∑_m∈ℒ𝐆_l,m P_l^t


    P_l^b=p̃'̃_l-Δ p'_l+p_l^E+∑_m∈ℒ𝐆_l,mP_l^t



    DP(t)=∑_l∈ℒ(p̃'̃_l(t)+r_lf̃_l(t))





 
 
 However, we can intuitively expect that having eq:error expression2 cannot be compliant with the mathematical formulation of the power flow equations (i.e., eq:complex power balance up-eq:relaxed internal current) in the case of insufficient capacity of ESS to satisfy the LDER constraint eq:error reduction rate. 
 
  The insufficient power rating of ESS means that possible p_l^E is small, making ϵ_l too large to comply with LDER constraint (see eq:error expression2). However, instead of making the problem infeasible, LDER constraint and eq:error expression2 are both satisfied by reducing the total dispatch error (left-hand side of eq:error expression2). This leads to the violation of the physical law of power flow because the grid losses error should take an unrealistic value that has the same order of magnitude as the prosumption deviation. In this way, the prosumption error and grid losses error cancel out each other to make the overall value of the left-hand side as small as the right-hand side of eq:error expression2. It would induce the increase of the current f_l such that the left-hand side of eq:relaxed internal current becomes strictly greater than the right-hand side, which leads to the inexactness of the solution. Therefore, we introduce an iterative algorithm, algo:losses realization, comprising two additional slack variables, γ_ϕ t^m and ζ_ϕ t, to replace eq:error expression2 by eq:modified error expression1 and eq:modified error expression2 such that the value of the internal current would never deviate away from the real value (i.e., the exactness of the solution is guaranteed). 
  	γ_ϕ t^m represents the realized grid losses deviation at mth iteration for scenario ϕ and time t, where m∈ℳ is the index of iterations of the algorithm. ζ_ϕ t indicates the unrealized part of the grid losses that should be updated to γ_ϕ t^m after each iteration for scenario ϕ and time t. γ_ϕ t^m achieves the accurate value of grid losses deviation, as the absolute value of ζ_ϕ t  reaches value below the defined tolerance (i.e., 1e-5 pu.).

    ∑_l∈ℒ(Δ p'_l ϕ t+r_lΔ f_l ϕ t)=∑_l∈ℒ(ϵ_l ϕ t+p^E_l ϕ t)+ζ_ϕ t,  ∀ϕ, ∀ t
    
       ∑_l∈ℒΔ p'_l ϕ t+γ_ϕ t^m =∑_l∈ℒ(ϵ_l ϕ t+p^E_l ϕ t),  ∀ϕ, ∀ t


Finally, the subproblem is described with an objective of minimization of the total grid losses and unserved load to satisfy the LDER constraint, and the operation period spans all days grouped into each day-type over the planning horizon.


    _∀φ',η,UL :   𝒮𝒞_d = Y N_d ∑_t ∈𝒯∑_ϕ∈Φ_dλ_ϕ d(w_l∑_l∈ℒr_lf_l ϕ t
    
    +w_u∑_l∈ℒ(ulp_l ϕ t^++ulp_l ϕ t^-+ulq_l ϕ t^++ulq_l ϕ t^-))



subject to: eq:error reduction rate, eq:active prosumption with unserved load, eq:reactive prosumption with unserved load eq:dispatch plan2, eq:AROPF_ESS_constraints-eq:modified error expression2,

    R_l = R_l^* : μ_ld,   C_l = C_l^* : ϑ_ld, ∀ l.
    
       Γ_d^(n)=[𝒮𝒞_d^*-∑_l∈ℒ(μ_l d(R_l-R_l^*)-ϑ_l d(C_l-C_l^*))], ∀ d, ∀ n

where UL={ulp^+,ulp^-, ulq^+,ulq^-} is the set of variables related to the unserved load. w_l and w_u are the weight coefficients associated with the grid losses minimization and unserved load, respectively. (<ref>) describes that the ESS power ratings and the energy reservoirs are fixed to the optimal solution values of the master problem. μ_l,d and ϑ_l,d are the duals of constraints related to the fixed ESS capacities, and are used to form the benders cuts for the master problems as shown in eq:bender's_cut. UB is calculated summing the optimal investment cost and the subproblem costs (i.e., UB=IC^*+∑_d∈𝒟𝒮𝒞_d^*).


§ SIMULATION




We validate the performance of the proposed methods with an existing Swiss distribution network with 55 bus and large capacity of renewable generation, as shown in fig:Aigle grid. The base voltage is 21kV and the base 3 phase power is 6MVA. 2.7MWp of PV generation capacity and 805kVA of hydropower generation capacity is installed. The planning horizon is set as 10 years, and we assume that the load consumption does not grow over the planning horizon. According to the indications of the operator of this grid, the number of candidate nodes for ESS installation is set as 5 out of 55 nodes (see tb:Simulation Parameters). To make sure to achieve the dispatchability, the penalty cost for the dispatch error is assumed as $700/MWh, which is significantly higher than a typical price settled in energy markets. The optimization problems are solved using the solver MOSEK via the MATLAB interface YALMIP. 


















	








































We have considered 8 typical day-types to cover the seasonal variation of the prosumption over the year. For each day-type, we assumed that prosumption forecasts for the simulation have been given from a reliable forecasting methodology. 1000 prosumption scenarios for each day are generated with equal probabilities based on the assumption that the prosumption profile follows a normal distribution as discussed in sec:achieving dispatchability. Then, we applied a K-medoids clustering <cit.> based scenario reduction technique to mitigate the computational burden. The number of reduced number of scenarios is determined by an algorithm explained in appendix:scenario. 


 §.§ Planning with 1 day under hourly dispatch

To illustrate the role of the 1st and 2nd block of the problem, we demonstrate the result of the simplified simulation considering an hourly dispatch for 1 representative day in two cases: case 1 with the original ampacity of the lines specified for the given grid and case 2 with the ampacity of line between node 4 and 20 reduced from the original value (see fig:Aigle grid). 32 representative scenarios obtained from the scenario reduction algorithm in appendix:scenario are used as operating scenarios. tb:installation_result1 indicates the optimal allocation result of ESS. The optimal result from the 1st block specifies the capacity in size of power rating and energy reservoir of ESS. However, the determined allocation of ESS cannot be guaranteed to be feasible and optimal to satisfy the LDER for the real operation of the grid,  since the grid losses and the ampacity constraint were neglected in the 1st block problem. After the 2nd block of the problem, as shown in the result of case 1 in tb:installation_result1, the optimal site of the ESS considering the objective of minimization of the grid losses is determined as node 4 and 27, resulting in the reduction of the grid losses and unserved energy compared to the result of the 1st block problem. 

In case 2, we can observe that a part of the load consumption was unserved to satisfy the LDER in the condition of restrained ampacity limit with the determined ESS allocation from the 1st block. In this regard, the result of the 2nd block shows the change in the allocation of the ESS due to the bottleneck of the line. The power rating and energy reservoir of the ESS on Node 4 is reduced, and those of the ESS on Node 27 is increased. tb:simplified_result shows that the unserved energy decreased to near zero after re-allocating the ESS.



	


     
	

	


 §.§ Planning with different number of scenarios
	




