Robust planning for distributed energy storage systems considering location marginal prices of distribution networks

Abstract

Energy storage plays an important role in integrating renewable energy sources and power systems, thus how to deploy growing distributed energy storage systems (DESSs) while meeting technical requirements of distribution networks is a challenging problem. This paper proposes an area-to-bus planning path with network constraints for DESSs under uncertainty. First, a distribution location marginal price (DLMP) formulation with maximum fluctuation boundaries of uncertainties is designed to select vulnerable areas exceeding voltage limits and higher line losses that occur in distribution networks. Different from simple multi-scenario power flow calculation and sensitivity analysis, DLMP with time and regional characteristics could be more intuitive to reflect line losses and voltage limits of distribution networks through price signals. After that, a two-stage stochastic robust optimization based planning method is developed to determine locations and capacities of DESSs in vulnerable areas. To make the uncertainty problem more tractable, stochastic scenarios are used to portray upper and lower boundaries of uncertainties, which avoids too-conservative decisions for robust optimization. Finally, numerical tests are implemented to testify the reasonability and validity of the proposed area-to-bus planning path under uncertainty. Compared with the DESSs planning framework without DLMP, the costs of DESSs are observably reduced with DLMP. With same budgets of uncertainty, investment costs of DESSs for the stochastic robust optimization with 30 and 50 scenarios are 3.91% and 4.45% lower than classical adaptive robust optimization (ARO).

Keywords

Distributed energy storage systems (DESSs)distribution location marginal price (DLMP)stochastic robust optimization system planning

Introduction

Background and significance

Energy storage systems are seen as an important part of efforts to boost intermittent renewable energy consumption while ensuring the stable operation of energy systems. In recent years, various centralized energy storage stations have been massively built around the world, such as 250 MW gateway energy storage project in California, and 100 MW energy storage demonstration project in Shanghai and so on.¹ However, these centralized power stations not only cover a huge area, but also require high operation and maintenance costs, in particular they are clumsy in preventing the expansion of battery accidents. In contrast, distributed energy storage systems (DESSs) have gradually emerged as the focus of markets due to their small capacity and lower expenses. Through absorbing or releasing excess power quickly, DESSs increase the self-consumption of renewable energy onsite and can offer significant cost savings to the private owner in distribution networks.^2,3 Most existing DESSs are not directly controlled by distribution system operators (DSOs), which could increase electricity prices at peak demand, influencing consumers.⁴ In other words, the planning for DESSs does not globally consider operating conditions of distribution networks, thus it is difficult to improve the power quality in vulnerable areas of distribution networks. Moreover, some repetitive and unreasonable planning for DESSs will not only increase investment costs but also affect the stability of distribution networks.^5,6 Some planning procedures without considering uncertain factors may cause deployments for DESSs divorced from reality. This paper attempts to answer the following key questions: 1) how to configure DESSs reasonably and orderly in vulnerable areas of distribution networks? 2) how to design the locations and capacities for DESSs scientifically under uncertainty?

Literature review

In order to enhance the flexibility of distribution networks in higher penetration of renewable energy sources, DESSs planning mostly revolves around load management,⁷ mitigation of voltage deviation,^8,9 peak-load shaving^10,11 and so forth. Researchers⁷ ascertain the optimal planning framework for battery energy storage to minimize network losses in terms of different load management. A scale of DESSs placement (e.g. uniform and non-uniform energy storage systems sizes) is developed to reduce voltage deviations and line losses.⁸ In order to delay distribution networks upgrade and mitigate voltage deviations, residual values, investment and operating costs of DESSs are integrated into a planning model.⁹ In Pimm et al.,¹⁰ the co-location of energy storage systems with loads could provide balancing services to grid, in which the potential of electricity storage for peak shaving is fully exploited. Similarly, an optimal configuration approach for battery units is applied to alleviate power unbalance, while achieving the maximum allowable photovoltaic capacity.¹¹ However, some grid and market factors (e.g. power flow and electricity prices) are not taken into account in the above approaches, which results in most energy storage systems being placed in unreasonable areas. It is difficult to capture potential economic values of DESSs and meet the fundamental requirements of distribution networks.

Some advanced approaches have been studied for planning distributed energy sources considering grid and market factors. In a nodal transmission-constrained energy market, price makers and geographically dispersed battery systems are coordinated using locational marginal price (LMP) in a nonlinear optimization model.¹² It is observed that such operation in the transmission line is suitable, but not always, desirable in distribution networks. In Singh and Goswami,¹³ a node pricing allocation scheme is developed for planning location and capacity of distributed energy sources. However, directly adapting the transmission level LMP to distribution systems and ignoring high R/X ratio could lead to large errors in marginal prices due to the discrepancy between distribution systems and transmission systems. For this reason, distribution location marginal price (DLMP) is applied to plan distributed energy sources in distribution networks. DSOs determine DLMP of distribution systems by taking into account charging constraints to alleviate disorderly charging with high penetration of electric vehicles.^14,15 An extension of former DLMP model based on novel linearized optimal power flow is developed by Yuan et al.¹⁶ and Wang et al.,¹⁷ in which accurate price information can be obtained. Inspired by the aforementioned literature, this paper develops a DLMP-based location planning strategy for DESSs. First, a price-guiding DLMP model in distribution networks is established according to the upper and lower boundaries of robust intervals of uncertainties. Second, the resulting schemes guide the preliminary planning of DESSs locations to reduce line losses and relieve voltage deviations for buses.

After a preliminary selection of planning areas, the detailed planning locations and capacities of DESSs are key points. Many optimization-based planning approaches have been studied, taking account into random renewable energy and loads. Stochastic programing-based planning approaches are developed to select suitable locations and capacities for energy storage systems.¹⁸ Whereas, these approaches make a tremendous computational burden due to a large number of scenario simulations.^19,20 Alternatively, robust optimization uses uncertain boundaries to cover worst-case scenarios for risk aversion, which alleviates the calculation pressure caused by stochastic scenario production.²¹ In order to facilitate the high percentage of renewable energy, an optimal allocation scheme for energy storage systems using robust optimization is proposed. Simulation results show that robust optimization schemes could reduce expansion planning costs while delaying network reinforcements.²² Considering N-1 contingencies, an adaptive robust optimization (ARO) model is established to confirm an optimal investment location for expansion planning.²³ To obtain the planning solution with minimum investing and operating costs, an ARO-based planning approach for energy storage units and electrical vehicle charging stations in distribution networks is developed considering short-term and long-term uncertainties.²⁴ In order to avoid too-pessimistic decisions of ARO-based approaches, distribution robust optimization (DRO) combines the advantages of stochastic programing and robust optimization to derive an optimal solution under the worst uncertain probability distribution.²⁵ Yang et al.²⁶ describe the coordinated planning for energy storage units and grid connection lines using DRO to minimize investment costs under wind abandonment constraints, but ignoring power flow in distribution networks. To optimize the sizing of battery energy storage in the power system, Guo et al.²⁷ propose a dedicated power flow model, in which inexact probability distributions are encapsulated in Wasserstein-metric based ambiguous set for the random renewable generation. However, the Wasserstein-metric based ambiguous set for DRO cannot be solved directly and its solution efficiency is impacted by the size of specimens.^28,29 Combining stochastic programing and robust optimization techniques for practical planning, this paper employs stochastic scenarios to portray an uncertainty set to simplify problem-solving processes. Fluctuation boundaries of stochastic scenarios are used to build upper and lower boundaries of robust intervals, which represent dynamic deviations from nominal values.

Contributions

This paper aims at solving the DESSs planning problem for DSOs under uncertainty. An area-to-bus planning path for DESSs is developed to meet fundamental requirements of distribution networks while handling uncertainties of renewable energy sources and loads. Compared with traditional planning methods, the proposed method accurately reduces voltage deviations and line losses. The stochastic scenarios and the robust interval are used to portray uncertainties of the planning process, in which the proposed model allows the upper and lower fluctuation intervals to be different. Therefore, mainly contributions of this paper can be summarized as follows.

The area-to-bus planning path is designed to optimize the locations and capacities of DESSs for DSOs under uncertainty, which alleviates voltage deviations and reduces line losses in distribution networks.

At the area level, according to the maximum boundary of uncertainty, DLMP approach is used to find vulnerable areas of the distribution network to guide the installation of DESSs at candidate buses. Fluctuation boundaries of stochastic scenarios are applied to depict upper and lower boundaries of uncertainty set.

At the bus level, the DESSs planning problem under uncertainty can be formulated as a two-stage stochastic robust optimization, which refrains too-conservative or over-optimistic decision schemes for stochastic programing-based planning and ARO-based planning approach. The effectiveness of the proposed model is verified in different stochastic scenarios.

The rest of this article is organized as follows. The problem statement of DESSs planning is given in Section 2. Section 3 introduces the power flow model and the stochastic-robust uncertainty model. Section 4 and Section 5 give a mathematical model of the proposed DESSs planning framework and solution methodology. Case studies and numerical results are discussed in Section 6 and conclusions are drawn in Section 7.

Problem statements

In this paper, the set of candidate buses for locating DESSs in distribution networks is defined as $M = {1, \dots, M}$ . $N = {1, \dots, N}$ denotes the set of buses in distribution networks and let $M \subseteq N$ . It is assumed that all energy storage units are invested by the same entity. In general, DESSs planning is mainly for local electricity services, thus it is subjected to economy, environment, and power constraints. In order to support the safe and reliable operation of distribution networks, DESSs are better placed in some vulnerable areas (e.g. voltage magnitude at lower bound and higher line losses), not in only commercial and residential areas representing users’ willingness to install energy storage systems. Moreover, DESSs planning should be compatible with renewable energy generation and loads in any bus of distribution networks. Therefore, the configuration process of DESSs in distribution networks is essentially an area-to-bus planning including a preliminary area planning problem and a detailed locations and capacities planning problem under uncertainty. A simplified illustration for DESSs planning in distribution networks is shown like Figure 1.

Figure 1.

A simplified illustration for DESSs planning in distribution networks.

First, the preliminary area planning problem for DESSs could be described as a DLMP problem with maximum fluctuation boundary of uncertainties. Based on the upper and lower boundaries of robust intervals of uncertainties, the preliminary planning finds a vulnerable area, in which power sold/purchased to/from upper-level grids are considered as decision variables. Different from traditional modes of operation (e.g. voltage limits through equipment), power market is an effective means to manage and utilize the power in the distribution network. The application of marginal pricing concepts in the distribution network can also be seen as signal of operational costs.³⁰ The price signal of DLMP is used to DESSs planning to alleviate the possible voltage deviation in the future distribution network with high renewable energy sources.

Second, the detailed locations and capacities planning problem of DESSs for DSOs under uncertainty is a two-stage planning problem after the preliminary selection of planning areas. The first stage pursues optimal locations and capacities for DESSs at candidate buses; then the second stage designs suitable operation schemes for DESSs. Given the existence of uncertainties of renewable energy sources and loads, the two-stage planning problem can be described as a two-stage stochastic robust optimization. The decision of DESSs with or not at candidate buses is treated as binary variable, and continuing variables include DESSs capacity and operation decisions (e.g. energy capacity and power capacity of installed DESSs, charging and discharging power of DESSs).

Preliminaries

Power flow model

Energy storage systems are accessed to regional distribution networks and transmit their power through transmission lines, which will undoubtedly have an impact on directions of power flow in distribution networks. Thus, power flow constraints are crucial for the DESSs planning model. Consider distribution networks described by a tree $H (N, E)$ , in which $N$ represents a set of buses and $E$ represents a set of lines. For each bus i ( $i \in N$ , except slack bus $1 \in N$ ), child buses are denoted by $C_{i}$ and $H_{i}$ is a set of parent buses of bus i. Time horizon for distribution network operating normally is defined as 1 day and divided $| T | = 24$ into equal-size time slots, denoted by $T = {1, \dots, 24}$ . $V_{i, t}$ and $I_{ij, t}$ are used, separately, to denote the squared voltage magnitude at bus i ( $i \in N$ ) and the current magnitude on line ij ( $ij \in E$ ). The definition of line power flow can be expressed by:

S_{ij, t} = v_{i, t} i_{ij, t}^{*}, \forall i \in N, ij \in E, t \in T,

(1)

where all buses satisfy power balance:

\sum_{j \in C_{i}} (S_{ij, t} - | i_{ij, t} |^{2} z_{ij}) - \sum_{k \in H_{i}} S_{jk, t} + S_{j, t}^{in} = 0, \forall j \in N, t \in T .

(2)

According to $v_{i, t} - v_{j, t} = z_{ij} i_{ij . t}$ and $S_{ij, t} = P_{ij, t} + j Q_{ij, t}$ , equation (1) can be split the real and imaginary parts as follows:

P_{ij, t} = \frac{r_{ij} v_{i, t} (v_{i, t} - v_{j, t} \cos φ_{ij, t})}{r_{ij}^{2} + x_{ij}^{2}} + \frac{x_{ij} v_{i, t} v_{j, t} \sin φ_{ij, t}}{r_{ij}^{2} + x_{ij}^{2}},

(3)

Q_{ij, t} = - \frac{r_{ij} v_{i, t} v_{j, t} \sin φ_{ij, t}}{r_{ij}^{2} + x_{ij}^{2}} + \frac{x_{ij} v_{i, t} (v_{i, t} - v_{j, t} \cos φ_{ij, t})}{r_{ij}^{2} + x_{ij}^{2}} .

(4)

(1) Linearized optimal power flow for networks: To obtain vulnerable areas in distribution networks, it is inefficient or even impracticable to straightly solve high R/X ratio using a linear approach. Thus, the linearized optimal power flow model¹⁶ is employed. To develop the linearized power flow formulas, two reasonable hypotheses based on generalization conditions need to be created below:

Assumption 1: In distribution networks, the voltage drop is commonly very tiny between two adjacent buses. The voltage angle difference $φ_{ij, t}$ between buses i and j should hence be approximated sufficiently to zero that $\sin φ_{ij, t} \approx φ_{i, t} - φ_{j, t}$ , $\cos φ_{ij, t} \approx 1$ .

Assumption 2: In local distribution networks, the voltage magnitude should happen under a reasonable condition (i.e. 1.0 p.u.). Therefore, $| v_{i, t} | \approx 1.0$ p.u. and $| v_{j, t} | \approx 1.0$ p.u.

Based on the above hypotheses, equations (3) and (4) could be a step forward to simplifying, which denote the model of novel linearized power flow.

P_{ij, t} \approx \frac{r_{ij} x_{ij}}{r_{ij}^{2} + x_{ij}^{2}} \frac{v_{i, t} - v_{j, t}}{x_{ij}} + \frac{x_{ij}^{2}}{r_{ij}^{2} + x_{ij}^{2}} \frac{φ_{i, t} - φ_{j, t}}{x_{ij}},

(5)

Q_{ij, t} \approx - \frac{r_{ij} x_{ij}}{r_{ij}^{2} + x_{ij}^{2}} \frac{φ_{i, t} - φ_{j, t}}{x_{ij}} + \frac{x_{ij}^{2}}{r_{ij}^{2} + x_{ij}^{2}} \frac{v_{i, t} - v_{j, t}}{x_{ij}} .

(6)

Similarly, matrix form of power injection at all buses can be expressed as:

[\begin{matrix} P_{j, t}^{in} \\ Q_{j, t}^{in} \end{matrix}] = [\begin{matrix} B_{2}^{'} \\ - B_{1}^{'} \end{matrix} \begin{matrix} B_{1}^{'} \\ B_{2}^{'} \end{matrix}] [\begin{matrix} φ_{ij, t} \\ v_{i, t} \end{matrix}] = B_{S} [\begin{matrix} φ_{ij, t} \\ v_{i, t} \end{matrix}],

(7)

B_{1} (i, j) = \frac{r_{ij, t}}{r_{ij, t}^{2} + x_{ij, t}^{2}}, \forall i \in N, t \in T, i \neq j,

(8)

B_{2} (i, j) = \frac{x_{ij, t}}{r_{ij}^{2} + x_{ij, t}^{2}}, \forall i \in N, t \in T, i \neq j,

(9)

B_{1} (i, i) = \sum_{j \in N, j \neq i} \frac{r_{ij}}{r_{ij, t}^{2} + x_{ij, t}^{2}}, \forall i \in N, t \in T,

(10)

B_{2} (i, i) = \sum_{j \in N, j \neq i} \frac{x_{ij}}{r_{ij, t}^{2} + x_{ij, t}^{2}}, \forall i \in N, t \in T .

(11)

In equation (7), $B_{S}$ shows a constant matrix constituted by the system impedance. Equations (8)–(11) define two matrices $B_{1}$ and $B_{2}$ , in which $B_{1}^{'}$ and $B_{2}^{'}$ are the submatrices of $B_{1}$ and $B_{2}$ excluding first row and first column.¹⁷

(2) DistFlow model for networks: Linearized optimal power flow could satisfy the search of vulnerable areas, but it could cause capacity deviation at buses. In this paper, optimal power flow is applied in detailed planning for DESSs. A DistFlow model is formulated in network $H$ . In order to change equations (1) and (2) into a convex problem, the corresponding relaxed equations can be written as:

\sum_{j \in C_{i}} (P_{ij, t} - I_{ij, t} r_{ij}) - \sum_{k \in H_{i}} P_{jk, t} + P_{j, t}^{in} = 0, \forall j \in N, t \in T,

(12)

\sum_{j \in C_{i}} (Q_{ij, t} - I_{ij, t} x_{ij}) - \sum_{k \in H_{i}} Q_{jk, t} + Q_{j, t}^{in} = 0, \forall j \in N, t \in T,

(13)

V_{j, t} = V_{i, t} - 2 (r_{ij} P_{ij, t} + x_{ij} Q_{ij, t}) + (r_{ij}^{2} + x_{ij}^{2}) I_{ij, t}, \forall ij \in E, t \in T,

(14)

I_{ij, t} = \frac{P_{ij, t}^{2} + Q_{ij, t}^{2}}{V_{i, t}}, \forall i \in N, ij \in E, t \in T,

(15)

v_{\min}^{2} \leq V_{i, t} \leq v_{\max}^{2}, \forall i \in N, t \in T,

(16)

0 \leq I_{ij, t} \leq i_{\max}^{2}, \forall ij \in E, t \in T .

(17)

Equations (12) and (13) describe active and reactive power conservations of distribution networks. At time slot t, voltage magnitude at bus i and current magnitude on line ij, respectively, can be calculated by equations (14) and (15). Equations (16) and (17) guarantee that the bus voltage magnitude and the line current magnitude would not beyond permitted ranges. However, the above model is non-convex due to the division operation among different variables in constraint (15). An inequality constraint (18) is adopted to relax the original model into a second-order cone programing³¹ that could be efficiently solved by commercial optimization solvers.

I_{ij, t} \geq \frac{P_{ij, t}^{2} + Q_{ij, t}^{2}}{V_{i, t}}, \forall i \in N, ij \in E, t \in T .

(18)

Stochastic-robust uncertainty model

In general, the deviation between the nominal value and the actual value of uncertain factors in the robust optimization is assumed to be fixed value at time slot t. Nevertheless, such descriptions for robust intervals could be regarded as too-pessimistic due to the fact that fluctuation ranges of uncertain factors are different in upper and lower boundaries. In this paper, fluctuation boundaries of stochastic scenarios are used to build dynamic deviations from nominal values at time slot t.

First, the stochastic scenario is used to characterize uncertainties. Here, wind speed, solar irradiance and loads are supposed to follow Weibull distribution, Beta distribution and Gaussian distribution³² respectively. Latin hypercube sampling method is adopted to generate approximately random samples $S \in Ω$ , through extracting a sample distributed evenly over a sample space. $Ω$ is number of stochastic scenarios. Compared with the classical Monte Carlo simulation method, it spreads the sample points more equally across all possible values by dividing each input distribution into intervals of equal probability.³³ Second, simultaneous backward reduction method is applied to select scenarios $S^{*} \in Ω^{*}$ , with the highest possibility of occurrence automatically. $S_{| Ω |, t}^{*}$ is typical scenarios selected through the above stochastic scenario method. $Ω^{*}$ denotes the set of typical scenarios. These methods are well-known in stochastic programing, please refer to Li and Xu³⁴ for interest. The upper and lower boundaries of uncertainty interval are the maximum and minimum values of all scenarios at time slot t.

{\bar{P}}_{t}^{u} = \max_{S^{*} \in Ω^{*}, t \in T} {S_{1, t}^{*}, S_{2, t}^{*}, \dots, S_{| Ω^{*} |, t}^{*}},

(19)

{\underline{P}}_{t}^{u} = \min_{S^{*} \in Ω^{*}, t \in T} {S_{1, t}^{*}, S_{2, t}^{*}, \dots, S_{| Ω^{*} |, t}^{*}} .

(20)

An occurrence probability accompanies each scenario, and the sum of all probabilities is equal to 1 (i.e. $\sum_{Ω^{*} \in Ω} ρ_{| Ω^{*} |} = 1$ ). According to occurrence probability of each scenario, nominal values of uncertain factors can be expressed as the expected value for all scenarios.

{\hat{P}}_{t}^{u} = \sum_{S^{*} \in Ω^{*}, t \in T} ρ_{| Ω^{*} |} S_{| Ω^{*} |, t}^{*} .

(21)

Based on nominal and boundary values obtained from the aforesaid scenarios, the robust set with budget of uncertainty is used to capture uncertainties. A straightforward structure of this uncertainty set could be, but is not restricted, as follows:

\begin{matrix} U = {\begin{matrix} {\tilde{P}}_{t}^{u} : \sum_{t \in T} \max {\frac{{\tilde{P}}_{t}^{u} - {\hat{P}}_{t}^{u}}{{\bar{P}}_{t}^{u} - {\hat{P}}_{t}^{u}}, \frac{{\tilde{P}}_{i, t}^{u} - {\hat{P}}_{t}^{u}}{{\hat{P}}_{t}^{u} - {\underline{P}}_{t}^{u}}} \leq Γ_{u}, \\ {\tilde{P}}_{t}^{u} \in [{\underline{P}}_{t}^{u}, {\bar{P}}_{t}^{u}], t \in T \end{matrix}} \end{matrix}

(22)

where ${\hat{P}}_{t}^{u}$ is the nominal value of uncertain factors at time t. ${\bar{P}}_{t}^{u} - {\hat{P}}_{t}^{u}$ and ${\hat{P}}_{t}^{u} - {\underline{P}}_{t}^{u}$ represent the deviation between the nominal and boundary value for uncertain factors, and the actual value of uncertain factors ${\tilde{P}}_{t}^{u}$ could be varied allodially in uncertainty interval $[{\underline{P}}_{t}^{u}, {\bar{P}}_{t}^{u}]$ . $Γ_{u}$ refers to the budget for uncertain factors taking an integer value between 0 and $| T |$ . Its value confines the total number of time intervals when uncertain factors deviate from its nominal value, thus it is employed to regulate the conservativeness of robust optimization problems. In equation (22), the range of robust intervals can be obtained by generating scenarios. The budget of uncertainty is the sum of the maximum value between upper and lower boundaries.

Preliminary area planning mathematical formulations

In this section, vulnerable areas for DESSs are selected based on the DLMP approach. At present, the safety of distribution networks is mainly considered congestions and line losses, however, the voltage limit is an urgent problem for radial distribution networks. First, the linearized optimal power flow model¹⁶ is employed to minimize operation costs of distribution networks. Then, a price-guiding DLMP formulation can be derived from the resulting power flow. DLMP with time and regional characteristics could be more intuitive to reflect line losses and voltage limits of distribution networks through price signals. Based on the price information, vulnerable areas can be determined to guide the installation of DESSs at candidate buses. The detailed mathematical formulation is listed as follows:

Preliminary area planning model

(A) : \min F^{dn} (z) = c_{a}^{b} P_{t}^{b} + c_{r}^{b} Q_{t}^{b} + c_{a}^{s} P_{t}^{s} + c_{r}^{s} Q_{t}^{s},

(23)

s . t . P_{t}^{b} + \sum_{i \in N} {\tilde{P}}_{i, t}^{res} = P_{t}^{s} + \sum_{i \in N} {\tilde{P}}_{i, t}^{l} + P_{t}^{loss}, \forall i \in N, t \in T,

(24)

Q_{t}^{b} = \sum_{i \in N} Q_{i, t}^{l} + Q_{t}^{s} + Q_{t}^{loss}, \forall i \in N, t \in T,

(25)

v_{\min} \leq v_{i, t} \leq v_{\max}, \forall i \in N, t \in T,

(26)

0 \leq P_{t}^{b} \leq P_{\max}^{b}, \forall t \in T,

(27)

0 \leq P_{t}^{s} \leq P_{\max}^{s}, \forall t \in T,

(28)

where $P_{t}^{loss} / Q_{t}^{loss}$ denote the total active/reactive power loss at time slot t, respectively. $z = {P_{t}^{b}, Q_{t}^{b}, P_{t}^{s}, Q_{t}^{s}}$ is the set of decision variables. Operation costs consist of the cost of power sold/purchased to/from upper-level grids. Equations (24) and (25) describe active and reactive power balance constraints of bus i, respectively. Equation (26) guarantees that the bus voltage would not violate allowable ranges. Equations (27) and (28) restrict the energy exchange between local distribution networks and upper-level grid.

$P^{loss}$ and $Q^{loss}$ can be expressed as a function of $P_{ij, t}^{2}$ , $Q_{ij, t}^{2}$ and $v_{i, t}^{2}$ , thus, not linearly associated with ${\tilde{P}}_{i, t}^{res}$ and ${\tilde{P}}_{i, t}^{l}$ . The linearization of $P^{loss}$ and $Q^{loss}$ could also be approximated by the first-order Taylor expansion, as expressed in Equations (29) and (30).

\begin{matrix} P_{t}^{loss} \approx P_{t}^{loss} ({\tilde{P}}_{i, t}^{res *}, {\tilde{P}}_{i, t}^{l *}) + \sum_{i \in N} \frac{\partial P_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{res}} ({\tilde{P}}_{i, t}^{res} - {\tilde{P}}_{i, t}^{res *}) \\ + \sum_{i \in N} \frac{\partial P_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{l}} ({\tilde{P}}_{i, t}^{l} - {\tilde{P}}_{i, t}^{l *}), \forall i \in N, t \in T, \end{matrix}

(29)

\begin{matrix} Q_{t}^{loss} \approx Q_{t}^{loss} ({\tilde{P}}_{i, t}^{res *}, {\tilde{P}}_{i, t}^{l *}) + \sum_{i \in N} \frac{\partial Q_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{res}} ({\tilde{P}}_{i, t}^{res} - {\tilde{P}}_{i, t}^{res *}) \\ + \sum_{i \in N} \frac{\partial Q_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{l}} ({\tilde{P}}_{i, t}^{l} - {\tilde{P}}_{i, t}^{l *}), \forall i \in N, t \in T . \end{matrix}

(30)

In order to obtain the relationship between voltage magnitudes and marginal prices for buses, equation (31) provides a detailed derivation of the linear function $v_{i, t}$ , which could be resulted in the associated pricing component about DLMP.

\begin{array}{l} v_{i, t} = \sum_{\begin{matrix} j \in N, \\ j \neq 1 \end{matrix}} [B_{N} (N B + i - 2, j - 1) P_{j, t}^{i n}] \\ + \sum_{\begin{matrix} j \in N, \\ j \neq 1 \end{matrix}} [B_{N} (N B + i - 2, N B + j - 2) Q_{j, t}^{i n}] \\ + \sum_{j \in P} \frac{v_{1, t}}{r_{1 j}^{2} + x_{1 j}^{2}} [B_{N} (N B + i - 2, j - 1) \\ + B_{N} (N B + i - 2, N B + j - 2)] \end{array}

(31)

where $v_{1, t}$ refers to the voltage magnitude of bus 1 and $j \in P$ denotes the bus j connected to bus 1. $x_{1 j} / r_{1 j}$ represent the reactance/resistance of the line connecting bus 1 and bus j, respectively. $B_{N}$ is the inverse matrix of $B_{S}^{- 1}$ , which could be expressed by equation (7).

Selecting vulnerable areas based on DLMP

Operation costs for distribution networks are not separable because of the existence of strong coupling relationship between constraints (24), (25), (26), (27) and (28), therefore which can be added into objective function using Lagrangian relaxation under uncertainty. The Lagrange function can be defined as follows:

\begin{matrix} (A_{L}) : {\begin{matrix} \min L (z) = F^{dn} (z) + λ_{p} p (z) + λ_{q} q (z) + μ_{i} g_{i} (z), \\ s . t . (24), (25), (26), (27), (28), \end{matrix} \end{matrix}

(32)

where $λ_{p}$ , $λ_{q}$ and $μ_{i}$ are the Lagrangian multipliers for corresponding constraints.³⁵ $p (z)$ and $q (z)$ represent equality constraints (24) and (25). $g_{i} (z)$ denotes inequality constraints (26)−(28). $A_{L}$ is a convex issue with a unique optimum. Karush - Kuhn - Tucker (KKT) conditions provide the optimality for the original problem and give an economical explanation of Lagrange multipliers. Based on Wang et al.,¹⁷ Lagrangian multipliers are marginal prices for analyzing the cost of power quality problems (e.g. voltage magnitude and higher line losses) in distribution networks. Therefore, the detailed definition of DLMP could be described as follows in equation (33) deriving from the KKT conditions, through which DSOs could check whether the energy storage planning is helpful to improve operational benefits. It should be noted that the DLMP for reactive power is not discussed in this paper, assuming the reactive power could be adjusted locally by users.

\begin{matrix} DLMP = \underset{ECC}{\underset{︸}{λ_{p}}} + \underset{LCC}{\underset{︸}{λ_{p} (\frac{\partial P_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{l}} + \frac{\partial P_{t}^{loss}}{\partial {\tilde{P}}_{i, t}^{res}})}} \\ + \underset{VCC}{\underset{︸}{\sum_{j \in N, j \neq 1} (- μ_{v_{j}}^{upper} + μ_{v_{j}}^{lower}) B_{N} (N_{l} + j - 2, i - 1)}} . \end{matrix}

(33)

In equation (33), $λ_{p}$ represent shadow prices of both equality constraints (24). $μ_{v_{j}}^{upper}$ and $μ_{v_{j}}^{lower}$ are shadow prices of inequality constraints (26). $μ_{v_{j}}^{upper}$ and $μ_{v_{j}}^{lower}$ are nonzero when $v_{i}$ breaches upper and lower boundaries of voltage magnitudes. The DLMP consists of three elements including energy cost component (ECC), voltage cost component (VCC) and loss cost component (LCC), which is related to shadow price $λ_{p}$ , $μ_{v_{j}}^{upper}$ and $μ_{v_{j}}^{lower}$ . DLMP with time and regional characteristics plays a significant role for helping DSOs identify vulnerable areas (e.g. marginal price over $ $α$ /MWh).

The DLMP is used to estimate the cost of exceeding voltage limits and higher line losses that occur in regional distribution networks. The marginal prices usually run within a reasonable range, which is given by DSOs based on historical experiences. Therefore, vulnerable areas are formed by areas in which the DLMP exceeds a reasonable price range. In fact, the apparent extent of vulnerable areas in distribution networks is subject to the worst-case uncertainties of renewable energy sources and loads. In constraint (22), as $Γ_{u}$ increases, uncertain factors could have more worst cases during scheduling horizon (i.e. $Γ_{u}$ confines the number of the worst-case realization of uncertainties). If the budget of uncertainty is higher, then the DLMP fluctuates more sharply, thus vulnerable areas are more obvious and their boundaries are clearer. To make DESSs planning beneficial to distribution networks, we try to find the greatest vulnerable areas. Thus, we select $Γ_{u}$ =24, the corresponding worst cases for renewable energy sources and loads are boundaries of robust interval. ${\tilde{P}}_{t}^{res}$ has upper boundary ${\bar{P}}_{i, t}^{res}$ and lower boundary ${\underline{P}}_{i, t}^{res}$ , ${\tilde{P}}_{t}^{l}$ has upper boundary ${\bar{P}}_{i, t}^{l}$ and lower boundary ${\underline{P}}_{i, t}^{l}$ . We have four combinations for these boundaries, that is, ${{\underline{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}$ , ${{\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}$ , ${{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}$ and ${{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}$ . Each combination corresponds to a unique vulnerable area like $DLM P_{{{\underline{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} \Rightarrow Ψ_{{{\underline{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , $DLM P_{{{\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}} \Rightarrow Ψ_{{{\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ , $DLM P_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} \Rightarrow Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , $DLM P_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}} \Rightarrow Ψ_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ Using the union operation for all vulnerable areas, one of the biggest vulnerable areas can be obtained in the preliminary planning process:

Ψ_{DLMP} = Ψ_{{{\underline{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} \cup Ψ_{{{\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}} \cup Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} \cup Ψ_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}} .

(34)

Solution methodology for the vulnerable area

The linearized optimal power flow model is calculated with Matpower 7.0.³⁶ Pseudo codes of acquired procedures for vulnerable areas using DLMP are outlined as Algorithm 1. In general, Algorithm 1 can be divided into three parts: Line 1 to line 2 is calculated DLMP in distribution networks using the linearized optimal power flow. Line 3 to line 6 found the corresponding vulnerable areas in four cases. Line 7 uses the union operation to obtain the biggest vulnerable areas.

Algorithm 1 Acquire vulnerable areas using DLMP
Input: $P_{t}^{b}, P_{t}^{s}, Q_{t}^{b}, Q_{t}^{s}, v_{i, t}, {\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}, {\underline{P}}_{i, t}^{res}, α$ .
1: Substitute four combinations for boundaries into model (23)−(28) to calculate DLMP for all buses,
2: ${DLMP}_{{{\underline{P}}_{i, t}^{r e s}, {\underline{P}}_{i, t}^{l}}}$ , ${DLMP}_{{{\underline{P}}_{i, t}^{r e s}, {\bar{P}}_{i, t}^{l}}}$ , $DLM P_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , $DLM P_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ .
3: if $DLMP > α$ then
4: Find the corresponding vulnerable areas in four cases,
5: $Ψ_{{{\underline{P}}_{i, t}^{r e s}, {\underline{P}}_{i, t}^{l}}}$ , $Ψ_{{{\underline{P}}_{i, t}^{r e s}, {\bar{P}}_{i, t}^{l}}}$ , $Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , $Ψ_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ .
6: end if
7: Use the union operation (34) to obtain the biggest vulnerable areas $Ψ_{DLMP}$ .
Output: Vulnerable areas to detailed planning of DESSs.

Detailed planning mathematical formulations

Objective function of DESSs planning

After preliminary planning to select vulnerable areas, the detailed planning locations and capacities of DESSs at candidate buses in vulnerable areas are key points. Given uncertainties of renewable energy sources and loads, this paper designs a two-stage stochastic robust optimization based DESSs planning approach to pursue the optimal DESSs planning scheme, in which voltage deviations for buses and line losses are fully considered. All buses in the vulnerable area have the potential to be candidate buses, thus, we consider that $M$ is derived from the vulnerable area $Ψ_{DLMP}$ . Detailed planning mathematical formulations for DESSs are expressed in the following.

\begin{matrix} (B) : F^{dp} (x, y, u) = \min_{x_{i}^{dess}, E_{i}^{dess}, P_{i}^{dess}} C^{inv} + \max_{{\tilde{P}}_{i, t}^{res}, {\tilde{P}}_{i, t}^{l}} \min_{\begin{matrix} P_{t}^{b}, P_{t}^{s}, P_{i, t}^{ch}, \\ P_{i, t}^{dis}, I_{ij, t}, V_{i, t} \end{matrix}} C^{om} \\ + C^{grid} + C^{loss} + D^{dev}, \end{matrix}

(35)

s . t . C^{inv} = \frac{σ {(1 + σ)}^{τ}}{{(1 + σ)}^{τ} - 1} \sum_{i \in M} (c_{E}^{inv} E_{i}^{dess} + c_{P}^{inv} P_{i}^{dess}),

(36)

C^{om} = c^{ch} \sum_{t \in T} η^{ch} P_{i, t}^{ch} - c^{dis} \sum_{t \in T} \frac{1}{η^{dis}} P_{i, t}^{dis},

(37)

C^{gd} = c_{a}^{b} \sum_{t \in T} P_{t}^{b} - c_{a}^{s} \sum_{t \in T} P_{t}^{s},

(38)

C^{ls} = c^{ls} \sum_{t \in T} \sum_{ij \in E} I_{ij, t} r_{ij},

(39)

D^{dev} = c^{ls} dw \sum_{i \in N} \sum_{t \in T} | \frac{v_{i, t} - v_{N}}{v_{N}} |,

(40)

0 \leq E_{i}^{dess} \leq x_{i}^{dess} E_{\max}^{dess}, \forall i \in M,

(41)

0 \leq P_{i}^{dess} \leq x_{i}^{dess} P_{\max}^{dess}, \forall i \in M,

(42)

\sum x_{i}^{dess} \leq n_{dess}, \forall i \in M,

(43)

0 \leq P_{i, t}^{ch} \leq P_{i, \max}^{ch}, \forall i \in M, t \in T,

(44)

0 \leq P_{i, t}^{dis} \leq P_{i, \max}^{dis}, \forall i \in M, t \in T,

(45)

SO C_{\min} \leq SO C_{i, t} \leq SO C_{\max}, \forall i \in M, t \in T,

(46)

S O C_{i, 1} = S O C_{i, 24} = 30 % S O C_{m a x}, \forall i \in ℳ,

(47)

SO C_{i, t} = SO C_{i, t - 1} + (η^{ch} P_{i, t}^{ch} - \frac{1}{η^{dis}} P_{i, t}^{dis}) E_{i}^{dess - 1}, \forall i \in M, t \in T,

(48)

P_{i, t}^{in} + P_{i, t}^{ch} + {\tilde{P}}_{i, t}^{l} + P_{t}^{s} = P_{i, t}^{dis} + {\tilde{P}}_{i, t}^{res} + P_{t}^{b}, \forall i \in N, t \in T .

(49)

$Eqs . (12) - (18), (27) - (28)$

where $x = {x_{i}^{dess}, E_{i}^{dess}, P_{i}^{dess}}$ refers to the first stage decision variables, in which $x_{i}^{dess}$ is binary variables for DESSs with or not at candidate buses. $y = {P_{i, t}^{ch}, P_{i, t}^{dis}, P_{t}^{b}, P_{t}^{s}, V_{i, t}, I_{ij, t}}$ is the second stage decision variables. $u$ denotes diverse uncertain factors, acquired by the uncertainty set $U$ (i.e. ${\tilde{P}}_{i, t}^{res} \in U_{res}$ , ${\tilde{P}}_{i, t}^{l} \in U_{l}$ ). In equation (35), the first stage is to minimize investment costs for all DESSs in distribution networks, and the second stage represents operation costs considering voltage deviations for buses and line losses. Equation (36) is annualized investment costs for DESSs. $σ$ denotes the annual discount rate and $τ$ represents the lifetime of DESSs, for example, 10 years. Equation (37) denotes operation costs for DESSs. Costs for the interactive power between distribution networks and upper-level grid could be expressed by equation (38). Equation (39) refers to the cost for line losses and the voltage deviation for bus is defined by equation (40), which can better reflect the improvement of power quality in the distribution network than voltage limits. d denotes the number of days in a year. w represents the conversion coefficient between the line loss and the voltage deviation for bus,³⁷ that is, $w = 500 MW / p . u .$ . $D^{dev}$ is converted into a price dimension in years using the unified conversion coefficient. Equations (41) and (42) restrict the energy and power for energy storage systems. equation (43) is the maximum installation number for energy storage systems, in which $x_{i}^{dess}$ is binary variable for DESSs with or not at the candidate buses. (44) and (45) limit the maximum charging and discharging power within the allowed range for DESSs, (46) enforces the state of charge at time slot t, (47) secures the identical dispatch maneuverability during each day and (48) represents the relationship between the state of charge and power for DESSs. In equation (49), the left-hand is the total power consumption for distribution networks, and the right-hand describes total available power generation for distribution networks, respectively.

Solution methodology

Currently, with cutting plane solution methods, like Benders decomposition method or column-and-constraint generation (C&CG) algorithm,³⁸ two stage robust optimization problems could be solved. Thus, C&CG algorithm is used to settle the planning model in the paper. First, the initial two stage DESSs planning model is transformed into a master problem (MP) determining DESSs with or not by discrete variable and a subproblem (SP) optimizing continuous system operation. MP and SP are then solved iteratively with discerning remarkable scenarios of uncertain factors in SP, and adding homologous constraints into MP. Eventually, MP and SP converge to an optimal solution by few iteratively solving. For a succinct exposition in this section, a compact matrix formulation for detailed planning of DESSs at candidate buses is given:

\min_{x} c^{T} x + \max_{u \in U} \min_{y \in Ω (x, u)} b^{T} y,

(50)

s . t . Ax \leq f, x \in {0, 1},

(51)

where

Ω (x, u) = {y : Dy \leq r,

(52)

Ey + Wu = g,

(53)

Fx + Gy \leq h,

(54)

{‖ K_{i} y + q ‖}_{2} \leq k_{i}^{T} y + p}

(55)

and x denotes investment decision variables including binary variables of the detailed planning at candidate buses. y refers to the continuous operation variables including DESSs operation decisions and other dependent variables. u represents uncertainties of renewable energy sources and loads captured by (22). Equation (43) can be converted to (51). Inequality constraints (16), (17), (27), (28), (44), (45), and (46) can be represented into (52). Equation (53) outlines the power balance that is, constraints (12), (13), (14), (15), (48), and (49). Equations (41) and (42) can be categorized as (54) that denotes the constraint with binary and continuous variables. Constraints (18) can also be written as (55).

(1) Master problem: MP of detailed planning at candidate buses is mostly developed to seek out the optimal investment decision under the uncertain scenarios at each iteration. It also provides a lower bound for the premier model (50). MP is expressed as:

MP : \min_{x, η} c^{T} x + η,

(56)

s . t . Ax \leq f, x \in {0, 1},

(57)

η \geq b^{T} y_{l}, \forall l \leq L,

(58)

D y_{l} \leq r, \forall l \leq L,

(59)

E y_{l} + {Wu}_{l}^{*} = g, \forall l \leq L,

(60)

Fx + G y_{l} \leq h, \forall l \leq L,

(61)

‖ K y_{l} + q ‖_{2} \leq k^{T} y_{l} + p, \forall l \leq L,

(62)

where L represents the iteration number and $l = 1, \dots, L$ . An auxiliary variable $η$ in MP is introduced. In MP of detailed planning, $y_{l}$ represents operating decisions, as new recourse variables at one per iteration of the algorithm appended to MP. $u_{l}^{*}$ denotes the optimal value of the uncertain parameter procured from SP at everytime iteration, hence, it is deemed fixed in MP. Note that MP of detailed planning is a mixed integer linear program, which can be immediately solved by existing standard branch-and-cut method.

(2) Subproblem under uncertainty: The goal of SP is to generate the worst-case scenarios of uncertainties for the aforementioned MP and likewise provide an upper bound on (50). SP is expressed as:

SP : \max_{u \in U} \min_{y \in Ω (x^{*}, u)} b^{T} y,

(63)

s . t . Dy \leq r,

(64)

Ey + Wu = g,

(65)

F x^{*} + Gy \leq h,

(66)

‖ Ky + q ‖_{2} \leq k^{T} y + p

(67)

where the investment decisions of DESSs at candidate buses are fixed as $x^{*}$ , which are iteratively acquired from the solution of MP at the algorithm. The worst-case scenarios of renewable energy sources and loads u and the homologous optimal DESSs operation y are found by solving the aforementioned subproblems. Based on the below KKT conditions, the max-min subproblem is therefore effectively converted to an equivalent max program.

KKT - SP : \max_{\begin{matrix} u \in U, \\ y \in Ω (x^{*}, u) \end{matrix}} b^{T} y,

(68)

s.t. (64)−(67)

D^{T} ϱ + E^{T} μ + G^{T} ω + K^{T} γ - k π = b,

(69)

(Dy - r)_{m} ϱ_{m} = 0,

(70)

(F x^{*} + Gy - h)_{n} ω_{n} = 0,

(71)

(‖ o ‖_{2} - ρ)_{i} κ_{i} + (o - Ky - q)_{i} γ_{i} + (ρ - k^{T} y - p)_{i} π_{i} = 0,

(72)

(D^{T} ϱ + E^{T} μ + G^{T} ω + K^{T} γ - k π - b)_{q} y_{q} = 0,

(73)

‖ γ ‖_{2} \leq π, κ = π,

(74)

ϱ, ω \leq 0, μ is free,

(75)

where $o = K_{i} y + q_{i}, ρ_{i} = k_{i}^{T} y + p_{i}$ . Equations (64)–(67) are primal constraints in detailed planning. (69) denotes dual constraints about (63)−(67) in SP of detailed planning. $ϱ$ , $ω$ , and $μ$ are the dual variables associated with linear constraints (64)−(66), respectively. $κ$ , $γ$ , and $π$ are dual variables of second-order cone constraint (67). Note that these complementary slackness conditions in constraints (70)−(73) are nonlinear, which can be linearized by employing the big-M method for bringing in binary variables. As an example in the case of (70), it then could be expressively reformulated as:

λ_{m} \geq - ϖ_{m} M, (Dy - r)_{m} \geq - (1 - ϖ_{m}) M, ϖ_{m} \in {0, 1},

(76)

where M denotes an adequately large constant and $ϖ_{m}$ represents a vector of binary variables. The constraint (76) implies the condition that, if $ϖ_{m} = 0$ then $ϱ_{m} = 0$ , $Dy \leq r$ , on the contrary, if $ϖ_{m} = 1$ then $ϱ_{m} \leq 0$ , $Dy = r$ .

Based on the aforementioned MP and SP, pseudo codes of procedures to solve the detailed planning problem using C&CG are outlined as Algorithm 2, in which line 2 means to solve MP, line 3 represents solving SP under uncertainty, and lines 4–8 determine whether the conditions are met. The solving flow chart of C&CG algorithm is shown in Figure 2.

Algorithm 2 C&CG algorithm for detailed planning of DESSs
Input: $M, P_{i, t}^{ch}, P_{i, t}^{dis}, P_{t}^{b}, P_{t}^{s}, V_{i, t}, I_{ij, t}$ .
1: Initialize $UB = \infty$ , $LB = - \infty$ , $L = 1$ , feasibility tolerance $ε$ .
2: Solve investment of DESSs under master problem to derive an optimal solution ( $x_{L + 1}^{}$ , $η_{L + 1}^{}$ , $y^{1 }$ ,…, $y^{L }$ ) and update $LB = c^{T} x_{L + 1}^{} + η_{L + 1}^{}$ .
3: Solve the optimal DESSs operation subproblem and update $UB = \min {UB, c^{T} x_{L + 1}^{} + b^{T} y_{L + 1}^{}}$ .
4: if $UB - LB \leq ε$ then
5: Return the optimal solution $x_{L + 1}^{*}$ and terminate.
6: else
7: Generate new variables $y^{L + 1}$ added constraints to MP and update L = L+1 for next iteration, return to step 2.
8: end if
Output: Locations and capacities of DESSs in distribution networks.

Figure 2.

The flow chart of C&CG algorithm.

Case study

In this section, different case studies are given to assess the performance of the proposed area-to-bus planning for DESSs in regional distribution networks. All simulations are performed on a laptop with an Intel Core i7 of 2.2 GHz and 16 GB of memory. Relevant problems are solved utilizing MATLAB (version 2019a) with YALMIP toolbox and GUROBI (version 9.5.0) solver.

In this paper, a modified IEEE 33-bus distribution network with four renewable energy sources at bus 4, 15, 21, 30 $\in N$ and loads locate at other buses are developed to verify the effectiveness of our proposed approach; see Figure 3. The regional distribution network is connected to upper-level power grid through the slack bus 1 $\in N$ . Without loss of generality, network parameters such as resistance and reactance of line are not modified in IEEE 33-bus distribution network. Part experimental data is provided by China electric power research institute. For example, time-of-use prices for selling/purchasing electricity to/from the upper-level grid are shown in Figure 4. In order to capture uncertainties, the number of generated scenarios of photovoltaic generations, wind generations and loads are 1000, respectively. The number of reduced scenarios of photovoltaic generations, wind generations and load are 970, respectively. Thirty representative stochastic scenarios for photovoltaic generations, wind generations and load in bus 13 are generated to portray robust intervals of uncertainties, shown in Figures 5 and 6. The DLMP of different combinations in the robust interval are shown in Figure 7.

Figure 3.

A modified IEEE 33-bus distribution network.

Figure 4.

Purchasing and selling electricity prices.

Figure 5.

Robust intervals of renewable energy sources based on stochastic scenarios: (a) the output power of wind generation in bus 4, (b) the output power of wind generation in bus 15, (c) the output power of photovoltaic generation in bus 21, and (d) the output power of photovoltaic generation in bus 21.

Figure 6.

Robust intervals of load in bus 13 based on stochastic scenarios.

Figure 7.

Different DLMP results of the modified IEEE 33-bus: (a) $DLM P_{{{\underline{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , (b) $DLM P_{{{\underline{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ , (c) $DLM P_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , and (d) $DLM P_{{{\bar{P}}_{i, t}^{res}, {\bar{P}}_{i, t}^{l}}}$ .

Comparison of DESSs planning with and without DLMP

This paper compares the results of DESSs planning methods with and without DLMP. Different from the DESSs planning method with DLMP for selecting vulnerable areas, the DESSs planning method without DLMP is an untargeted method. According to equation (34), if $DLMP > α$ , then the corresponding vulnerable areas can be selected based on robust intervals. Table 1 presents the candidate bus from different vulnerable areas under different $α$ . Case 1 denotes the union of all vulnerable areas $Ψ_{DLMP}$ , when $DLMP > α = $ 45 / MWh$ . Case 2 is similar to that, representing $DLMP > α = $ 50 / MWh$ . Case 3 refers to a typical vulnerable area $Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ , when $DLM P_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} > α = $ 45 / MWh$ . Thus, $Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ in case 3 is a subset of $Ψ_{DLMP}$ in case 1.

Table 1.

Candidate bus in vulnerable areas with different DLMP.

Buses
Cases	Candidate bus in vulnerable areas
Case 1	7, 8, 9, 10, 11, 12, 14, 16, 17, 26, 27, 28, 29, 31, 32, 33
Case 2	7, 8, 9, 10, 11, 12, 14, 16, 17, 31, 32, 33
Case 3	8, 9, 10, 11, 12, 14, 16, 17

In the detailed planning, let $Γ_{wt} = Γ_{load} = 12$ and $Γ_{pv} = 6$ . Comparison of power profiles under different planning methods with DLMP and without DLMP are shown in Figure 8. The power profiles of DESSs are different in the case with DLMP and without DLMP due to different installation locations. Comparison of planning schemes is depicted clearly in Table 2, in which the computing time of the planning with DLMP is less than that of the planning method without DLMP planning. Meanwhile, the costs of DESSs in case 1 and case 2 are less than the case without DLMP. Compared to the total installed capacity of the case without DLMP, the capacities of case 1 and case 2 are decreased by 3.84% and 0.44% respectively. The total capacity of case 3 is increased by 3.13% compared to that of the case without DLMP. This is because case 3 can only find a part of the vulnerable area, thus the selected candidate bus in $Ψ_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}}$ might not satisfy planning requirements for the distribution network. One can conclude that DLMP could help the planning for DESSs reduce overall planning costs and improve planning effectiveness.

Figure 8.

Comparison of power profiles under different planning methods: (a) $DLMP > $ 45 / MWh$ , (b) $DLM P_{{{\bar{P}}_{i, t}^{res}, {\underline{P}}_{i, t}^{l}}} > $ 45 / MWh$ , (c) $DLMP > $ 50 / MWh$ , and (d) without DLMP.

Table 2.

Comparison results with and without DLMP methods.

Planning method	Without DLMP	With DLMP
Planning method	Without DLMP	Case 1	Case 2	Case 3
Location	2, 17, 28	7, 17, 26	7, 17, 33	8, 9, 17
Time (sec)	2260	1158	1140	1122
Investment costs ($)	3054.87	2984.86	3047.82	3130.50
Capacity (MW)	0.2778 (2) 0.2703 (17) 0.2685 (28)	0.2564 (7) 0.2703 (17) 0.2597 (26)	0.2778 (7) 0.2685 (17) 0.2667 (33)	0.2778 (7) 0.2941 (9) 0.2703 (17)
Total capacity (MW)	0.8166	0.7864	0.8130	0.8422

To demonstrate the effects of DESSs deployed in the distribution network, DLMPs of case 1 with DESSs and without DESSs are depicted in Figure 9. Compared to the DLMP in Figure 9(a), the DLMP in Figure 9(b) is significantly reduced. Especially at 21:00, the DLMP is reduced by $66.73/MWh. This is because the voltage offset at buses could be adjusted flexibly and the line loss can be reduced by charging and discharging of DESSs. Select the time point at 20:00 as the research point, voltage deviations for buses in case 1 with DESSs and without DESSs are shown in Figure 10. Voltage deviations for buses in the case with DESSs are lower than those in the case without DESSs. In bus 29, the voltage deviation is greatly reduced by 26.94%, which satisfies the power supply quality requirements of the distribution network. The price signal of DLMP is used to DESSs planning to alleviate the possible voltage deviation in the future distribution network with high renewable energy sources.

Figure 9.

The DLMP with DESSs and without DESSs: (a) the DLMP without DESSs and (b) the DLMP with DESSs.

Figure 10.

Voltage deviation for buses with and without DESSs.

Comparison of DESSs planning under different budgets of uncertainty

As mentioned in preliminaries, the conservativeness of robust optimization problems can be dominated by adjusting $Γ$ , in which $Γ$ takes integer values. In this part, we test these characteristics of the devised planning schemes of case 1 under different budget realizations of uncertainty, namely $Γ = 0$ , $Γ = 6$ , $Γ = 12$ .

Figure 11 shows that comparison of power profiles under different budgets of uncertainty. With the increase of uncertainty, DESSs are usually charged when the output of renewable energy sources is high and discharged when the power is short with the output of renewable energy sources insufficiently. The distribution network is motivated by electricity prices selling power to upper-level grids to reduce operation costs. As $Γ$ increases, larger deviations from the nominal value of the output of renewable energy sources and loads are considered, which means the output power of renewable energy is likely to become less and loads are higher. Therefore, the distribution network has properly scheduled DESSs to cover these power deviations.

Figure 11.

Renewable energy source and load profiles under different budgets of uncertainty: (a) $Γ = 0$ , (b) $Γ = 6$ , and (c) $Γ = 12$ .

Table 3 lists comparison results under different budgets of uncertainties. Although the locations of DESS planning varies with $Γ$ , some installation locations are still similar. As $Γ$ increases, the capacities of DESSs are expanded to handle uncertainty, thus the corresponding investment costs also gradually increase. It is observed that robust optimization takes the worst-case scenarios into decision-making processes, which guarantees the robustness of DESSs planning schemes. In the same way, a higher budget of uncertainties increases in worse economics but better against the capability of risk (i.e. higher conservatism level of the model).

Table 3.

Comparison results under different budgets of uncertainty.

$Γ$	Location	Investment costs ($)	Capacity (MW)	Total capacity (MW)
0	7,17,32	2940.22	0.2564 (7) 0.2564 (17) 0.2532 (32)	0.7660
6	7,17,26	2973.25	0.2634 (7) 0.2667 (17) 0.2564 (26)	0.7865
12	7,17,26	3001.07	0.2614 (7) 0.2703 (17) 0.2597 (26)	0.7914

Comparison of stochastic robust optimization and classical ARO

The different comparisons among stochastic robust optimization and classical ARO are shown in this section. Note that parameters selected by all optimization approaches are identical but the uncertainty sets are slightly different due to different upper and lower boundaries of the corresponding robust interval. In classical ARO, the uncertainty sets usually run within a reasonable range, which is given by historical experiences, as shown in Figures 12 and 13. Fifty representative stochastic scenarios for photovoltaic generations, wind generations and loads are generated to portray robust intervals of uncertainties, see Figures 14 and 15. The deviation of the nominal value from the upper and lower boundaries is usually equal in the classical ARO. However, stochastic robust optimization uses stochastic scenarios to describe the interval boundaries, and the deviation between the nominal value and the upper/lower boundaries may be unequal. Compared with classical ARO, stochastic robust optimization can better reduce conservatism.

Figure 12.

Interval forecasts of renewable energy sources: (a) the output power of wind generation in bus 4, (b) the output power of wind generation in bus 15, (c) the output power of photovoltaic generation in bus 21, and (d) the output power of photovoltaic generation in bus 30.

Figure 13.

Interval forecasts of load in bus 13.

Figure 14.

Robust intervals of renewable energy sources based on 50 scenarios: (a) the output power of wind generation in bus 4, (b) the output power of wind generation in bus 15, (c) The output power of photovoltaic generation in bus 21, and (d) THE output power of photovoltaic generation in bus 30.

Figure 15.

Robust intervals of load in bus 13 based on 50 scenarios.

In this part, let $Γ_{wt} = Γ_{load} = 12$ and $Γ_{pv} = 6$ . Figures 8(a) and 16 show that comparison of power profiles under stochastic robust optimization with 30/50 scenarios and classical ARO. The power profiles of DESSs are different in stochastic robust optimization with 30/50 scenarios and classical ARO due to different boundaries of robust interval. Table 4 reports the comparison results utilizing stochastic robust optimization and classical ARO. It can be observed that investment costs and capacities are gradually decreasing compared with classical ARO. This is because the model for stochastic robust optimization is less conservative than classical ARO, and the large number of stochastic scenarios further reduces the conservativeness for robust optimization. With the same budgets of uncertainty, investment costs of DESSs for the stochastic robust model with 30 and 50 scenarios are 3.91% and 4.45% lower than the classical ARO, respectively. Using more stochastic scenarios can be closer to reality, better describe the dynamic deviation from nominal values of uncertain factors, and reduce the conservativeness of DESSs planning. Therefore, the stochastic robust optimization avoids too-pessimistic decisions and addresses diverse uncertainties.

Figure 16.

Comparison of power profiles under stochastic robust optimization and classical ARO: (a) stochastic robust optimization with 50 scenarios and (b) classical ARO.

Table 4.

Comparison results for stochastic robust optimization and classical ARO.

Planning method	Classical ARO	Stochastic robust optimization
Planning method	Classical ARO	30 scenarios	50 scenarios
Investment costs ($)	3106.30	2984.86	2967.99
Capacity (MW)	0.2703 (7) 0.2941 (17) 0.2685 (26)	0.2564 (7) 0.2703 (17) 0.2597 (26)	0.2548 (7) 0.2667 (17) 0.2584 (26)
Total capacity (MW)	0.8329	0.7864	0.7799

Conclusions

An area-to-bus planning path is designed to obtain acceptable vulnerable areas of distribution networks in preliminary area planning as well as determine locations and capacities of DESSs at candidate buses in vulnerable areas with detailed planning.

Different from the DESSs planning method with DLMP for selecting vulnerable areas, the DESSs planning method without DLMP is an untargeted method. Simultaneously, DLMP with time and regional characteristics could be more intuitive to reflect line losses and voltage limits of distribution networks through price signals. Under same conditions of uncertainty, the total capacity of DESSs in the method with DLMP is less than that in the method without DLMP. The total capacity savings are 3.84%, 0.44%, respectively. Thus, DLMP could help the DESSs planning reduce overall planning costs and improve planning effectiveness. If the method only considers a part of the vulnerable area of robust interval, it may not satisfy planning requirements for the distribution network. Hence, the biggest vulnerable area needs to be derived considering different boundaries of uncertainty. The stochastic robust optimization can describe the dynamic deviation from nominal values of uncertainties. With the same budgets of uncertainty, investment costs of DESSs for the stochastic robust model with 30 and 50 scenarios are 3.91% and 4.45% lower than the classical ARO, respectively. Thus, the stochastic robust approach avoids too-pessimistic decisions than classical ARO approach.

The proposed DESSs planning framework conforms to the network constraints of power systems. Thus, the theoretical framework and implementation methods for DESSs planning are developed to provide strong support for DSOs in practical applications. In order to make the model tractable, this paper adopts stochastic robust optimization. However, the accuracy of uncertainty description is still lacking. Therefore, future work will be studied DRO to enhance its precision.

Footnotes

Appendix

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Yue Sun

References

Khezri

Mahmoudi

Aki

Optimal planning of solar photovoltaic and battery storage systems for grid-connected residential sector: review, challenges and new perspectives. Renew Sustain Energ Rev 2022; 153: 111763.

Zakeri

Gissey

Dodds

, et al. Centralized vs. Distributed energy storage – benefits for residential users. Energy 2021; 236: 121443.

Vidhya

Balaji

Modelling, design and control of a light electric vehicle with hybrid energy storage system for Indian driving cycle. Meas Control 2019; 52: 1420–1433.

Yin

Chen

Liang

, et al. Robust transmission network expansion planning based on a data-driven uncertainty set considering spatio-temporal correlation. Sustain Energy Grids 2023; 33: 100965.

Das

Bass

Kothapalli

, et al. Overview of energy storage systems in distribution networks: placement, sizing, operation, and power quality. Renew Sustain Energ Rev 2018; 91: 1205–1230.

Sun

, et al. Deep learning based semi-supervised control for vertical security of maglev vehicle with guaranteed bounded airgap. IEEE T Intell Transp 2021; 22(7): 4431–4442.

Zhang

Ren

Dong

, et al. Optimal placement of battery energy storage in distribution networks considering conservation voltage reduction and stochastic load composition. IET Gener Transm Distrib 2017; 11: 3862–3870.

Das

Bass

Kothapalli

, et al. Optimal placement of distributed energy storage systems in distribution networks using artificial bee colony algorithm. Appl Energy 2018; 232: 212–228.

Das

Bass

Mahmoud

, et al. An optimal allocation and sizing strategy of distributed energy storage systems to improve performance of distribution networks. J Energy Storage 2019; 26: 100847.

10.

Pimm

Cockerill

Taylor

PG.

The potential for peak shaving on low voltage distribution networks using electricity storage. J Energy Storage 2018; 16: 231–242.

11.

Wang

Zhang

Dong

, et al. Improving hosting capacity of unbalanced distribution networks via robust allocation of battery energy storage systems. IEEE Trans Power Syst 2021; 36: 2174–2185.

12.

Mohsenian-Rad

Coordinated price-maker operation of large energy storage units in nodal energy markets. IEEE Trans Power Syst 2016; 31: 786–797.

13.

Singh

Goswami

SK.

Optimum allocation of distributed generations based on nodal pricing for profit, loss reduction, and voltage improvement including voltage rise issue. Int J Elec Power 2010; 32: 637–644.

14.

Liu

Oren

, et al. Distribution locational marginal pricing for optimal electric vehicle charging through chance constrained mixed-integer programming. IEEE Trans Smart Grid 2018; 9: 644–654.

15.

Oren

SS.

Distribution locational marginal pricing for optimal electric vehicle charging management. IEEE Trans Power Syst 2014; 29: 203–211.

16.

Yuan

Wei

, et al. Novel linearized power flow and linearized OPF models for active distribution networks with application in distribution LMP. IEEE Trans Smart Grid 2018; 9: 438–448.

17.

Wang

Tang

, et al. Three-phase DLMP model based on linearized power flow for distribution with application to DER benefit studies. Int J Elec Power 2021; 130: 106884.

18.

Zheng

Hill

Dong

ZY.

Multi-agent optimal allocation of energy storage systems in distribution systems. IEEE Trans Sustain Energy 2017; 8: 1715–1725.

19.

Wang

Zhang

, et al. Robust distributed optimization for energy dispatch of multi-stakeholder multiple microgrids under uncertainty. Appl Energy 2019; 255: 113845.

20.

Zhang

Wang

, et al. Robust optimization for energy transactions in multi-microgrids under uncertainty. Appl Energy 2018; 217: 346–360.

21.

Yan

Geng

Bie

, et al. Two-stage robust energy storage planning with probabilistic guarantees: A data-driven approach. Appl Energy 2022; 313: 118623.

22.

Lyons

Davison

, et al. Robust scheduling scheme for energy storage to facilitate high penetration of renewables. IEEE Trans Sustain Energy 2016; 7: 797–807.

23.

Pinto

Unsihuay-Vila

Tabarro

FH.

Coordinated operation and expansion planning for multiple microgrids and active distribution networks under uncertainties. Appl Energy 2021; 297: 117108.

24.

Baringo

Boffino

Oggioni

Robust expansion planning of a distribution system with electric vehicles, storage and renewable units. Appl Energy 2020; 265: 114679.

25.

Gao

Wang

Liu

, et al. Data-driven distributionally robust joint planning of distributed energy resources in active distribution network. IET Gener Transm Distrib 2020; 14: 1653–1662.

26.

Yang

Xei

Wei

, et al. Coordinated planning of storage unit in a remote wind farm and grid connection line: A distributionally robust optimization approach. In: IEEE Innovative Smart Grid Technologies-Asia (ISGT Asia), Chengdu, China, 21 May-24 May 2019, pp.424–428. New York: IEEE.

27.

Guo

Wei

Chen

, et al. Sizing energy storage to reduce renewable power curtailment considering network power flows: a distributionally robust optimisation approach. IET Renew Power Gener 2020; 14: 3273–3280.

28.

Guo

Baker

Dall’Anese

, et al. Data-based distributionally robust stochastic optimal power flow—Part I: Methodologies. IEEE Trans Power Syst 2019; 34: 1483–1492.

29.

Duan

Fang

Jiang

, et al. Distributionally robust chance-constrained approximate AC-OPF with Wasserstein metric. IEEE Trans Power Syst 2018; 33: 4924–4936.

30.

Heydt

Chowdhury

Crow

, et al. Pricing and control in the next generation power distribution system. IEEE Trans Smart Grid 2012; 3: 907–914.

31.

Wang

Sun

, et al. Robust estimation of reactive power for an active distribution system. IEEE Trans Power Syst 2019; 34: 3395–3407.

32.

Wang

Zhang

Song

, et al. Stochastic cooperative bidding strategy for multiple microgrids with peer-to-peer energy trading. IEEE Trans Ind Inform 2022; 18: 1447–1457.

33.

Risk-averse coordinated operation of a multi-energy microgrid considering voltage/var control and thermal flow: an adaptive stochastic approach. IEEE Trans Smart Grid 2021; 12: 3914–3927.

34.

Temporally-coordinated optimal operation of a multi-energy microgrid under diverse uncertainties. Appl Energy 2019; 240: 719–729.

35.

Wang

, et al. Optimal distributed transaction of multiple microgrids in grid-connected and islanded modes considering unit commitment scheme. Int J Elec Power 2021; 132: 107146.

36.

Zimmerman

Murillo-Sanchez

Thomas

RJ.

MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education. IEEE Trans Power Syst 2011; 26: 12–19.

37.

Zulpo

Leborgne

Bretas

AS.

Optimal siting and sizing of distributed generation through power losses and voltage deviation. In: 16th International Conference on Harmonics and Quality of Power (ICHQP), Bucharest, Romania, 25 May-28 May 2014, pp.871–875. New York: IEEE.

38.

Zeng

Zhao

Solving two-stage robust optimization problems using a column-and-constraint generation method. Oper Res Lett 2013; 41: 457–461.