Characterizing and Visualizing the Impact of Energy Storage on Renewable Energy Curtailment in Bulk Power Systems

: The uncertain natures of renewable energy lead to its underutilization; energy storage unit (ESU) is expected to be one of the most promising solutions to this issue. This paper evaluates the impact of ESUs on renewable energy curtailment. For any ﬁxed renewable power output, the evaluation model minimizes the total amount of curtailment and is formulated as a mixed integer linear program (MILP) with the complementarity constraints on the charging and discharging behaviors of ESUs; by treating the power and energy capacities of ESUs as parameters, the MILP is transformed into a multi-parametric MILP (mp-MILP), whose optimal value function (OVF) explicitly maps the parameters to the renewable energy curtailment. Further, given the inexactness of uncertainty’s probability distribution, a distributionally robust mp-MILP (DR-mp-MILP) is proposed that considers the worst distribution in a neighborhood of the empirical distribution built by the representative scenarios. The DR-mp-MILP has a max–min form and is reformed as a canonical mp-MILP by duality theory. The proposed method was validated on the modiﬁed IEEE nine-bus systems; the parameterized OVFs provide insightful suggestions on storage sizing.


Introduction
The consumption of fuel fossils has considerably soared in the past few decades due to the industrialization and economic development worldwide, resulting in the issues of energy shortage and environmental pollution. Under this background, the sector of renewable energy generation (REG) springs up; by 2019, the global installed capacity of wind and solar generation, two dominating technologies of REG, expanded to 651 [1] and 627 GW [2], respectively. Such a high penetration level of REG will cause difficulties to the power system operation: On the one hand, the inherent uncertainty of renewable energy tests the supply reliability and operating security of the power grid [3,4]; on the other hand, the imbalance between supply and demand, as well as the transmission congestion, underutilizes the renewable energy and induces a lot of energy curtailment [5]. With the aim to reduce the renewable energy curtailment, researchers pin their hopes on the energy storage unit (ESU), which can store the excessive renewable energy and discharge it when needed, thus cutting down the spillage.
The planning of ESUs in support of reducing the renewable energy curtailment has been widely investigated. Dui X. et al. proposes a two-stage method to determine the optimal capacity of battery energy storage to decrease wind power curtailment in grid-connected wind farms [6]; the first stage schedules the unit commitment of thermal multi-parametric programming, which is a systematic theory and technique that subdivides the parameter set into characteristic regions where the optimal values are given as explicit functions of the parameters [27]; hence, if we treat the power and energy capacities as parameters, such functions provide the decision-maker with insightful information about the impact of ESU on renewable energy curtailment.
Multi-parametric programming has been used in some related works. In Reference [28], the impact of ESU on the operational economy of distribution network is discussed; similarly, Reference [29] focuses on the economic influence of ESU in energy hub; however, the complementarity constraints on the charging and discharging behaviors of ESUs are not imposed by binary variables in these two papers, which is not applicable in this paper, as illustrated in Section 2.3. The evaluation model in Reference [30] follows the scheme of robust optimization.
This paper combines the distributionally robust optimization and multi-parametric programming, and proposes a distributionally robust multi-parametric mixed-integer linear program (DR-mp-MILP) to assess the impact of energy storages on renewable energy curtailment under distributional uncertainty. The DR-mp-MILP has a max-min form and is reformed as a canonical mp-MILP by duality theory; the results from the IEEE nine-bus system tests validate the effectiveness and performance of our proposed method, and provide insightful suggestions on storage sizing.
The rest of this paper is organized as follows: In Section 2, the fundamental models of the power system are first introduced, followed by an illustrative example about the necessity to use binary variables to avoid simultaneous charging and discharging. In Section 3, the multi-parametric formulations are proposed. Section 4 presents the solution to mp-MILP, and Section 5 shows the results of case studies. Discussions are given in Section 6.

Basic Models
This section introduces the fundamental models for the problem under study. The operating constraints of power system are first given, followed by the objective function to optimize the renewable energy curtailment. The power flows in the network are described by the renowned direct-current power flow model, which is loss-less and neglects the reactive power distribution. We make some clarifications in the revised manuscript. Specifically, an example is designed to illustrative that it is inevitable to use binary variables to impose complementarity constraints on the discharging and charging behaviors of an ESU.

Operating Constraints
Considering the multi-period operation of a power system with renewable plants and ESUs, the constraints are gathered as shown below: 0 ≤ p r kt ≤ p rm kt , ∆p rm kt = p rm kt − p r kt , ∀k ∈ W, ∀t ∈ T where (1) describes the system-wise power balance. Moreover, (2) says the active power flow in line l cannot exceed its transmission capacity; π li is the power transfer distribution factor from bus i to line l; p n it represents the net power injection at bus i in period t, which depends on the production of generator (p g it ), the power exchange with ESU (p sc it /p sd it ), the dispatched renewable power (p r it ), as well as the load (p d it ), if there is any. Moreover, (3) restricts the output and ramping rates of generators.
The charging and discharging behaviors of ESUs are constrained by (4), where P s j is the power capacity; when the binary indicator is b s jt = 1/0, the ESU j operates under the charging/discharging mode in period t, with a corresponding power lower than P s j . The state-of-charge (SoC) dynamics of ESU is expressed by (5), where E s j is the energy capacity, and α l j is a constant coefficient indicating the lower bound of SoC; e j0 is the initial level of SoC, and h is the duration of one time slot, i.e., one hour; besides, the terminal SoC is set to the initial state e j0 , implied by (6).
In (7), p rm kt is the available renewable power in renewable plant k in period t, which is uncertain and determined by the weather conditions; hence, the dispatched renewable power, p r kt , cannot exceed p rm kt . Meanwhile, the curtailed renewable power is defined as ∆p r kt .

Objective Function
For the purpose of reducing the renewable energy curtailment, the objective function is as follows:

Illustration: Necessary Charging-Discharging Complementarity Constraint
For a storage-concerned economic dispatch problem, the complementarity constraints on charging and discharging are, in most cases, naturally satisfied without using binary variables. However, in this paper, such binary variables are necessary because, in the renewable energy curtailment problem, simultaneous charging and discharging can act as a load, which consumes a certain amount of renewable energy and thus reduces the curtailment. Here, we give a simple example for illustration: The system consists of one bus, one ESU, and one wind turbine. Let T = 1, P s = 1 MW, E s = 2 MWh, e 0 = 2 MWh, η c = η d = 0.9, p d = 0.6MW, p rm = 1MW; for brevity, the bus and time subscripts are omitted.
If we use the binary variable as (4), simultaneous charging and discharging is not permitted; meanwhile, the SoC level has reached its ceiling. Consequently, 0.6 MWh of wind energy is consumed by the load, and the rest 0.4 MWh is curtailed.
However, if we do not use the binary variable, consider the following strategy: p sc = 1 MW, p sd = 0.81 MW; then, e 1 = e 0 + h p sc η c − p sd /η d = 2 MWh. According to the power balance condition in (1), the dispatched wind power p r = 0.79 MW, and thus the curtailed wind energy, is 0.21 MWh. There is no doubt such a strategy is feasible while the curtailment is lower than that in the case with binary variable. This example verifies that the binary complementarity indicator is necessary in this paper; otherwise, simultaneous charging and discharging will happen.

Multi-Parametric Formulations
This section proposes the deterministic and distributionally robust mp-MILP formulations.

Deterministic mp-MILP
To give the closed form of deterministic mp-MILP, we first define the following vectors: • Uncertain variable vector:ŵ = p rm kt , ∀k ∈ W, ∀t ∈ T T .
• Continuous variable vector:x is a column vector and encloses all the remaining continuous variables in (1)- (7).
Then, any equality in (1)-(7) can be replaced by two opposite inequalities; for example, (6) is equivalent to the following: By far, the closed form of deterministic mp-MILP can be written as follows: whereÂ,B,b,Ĉ,F, andĉ are coefficient matrices or vectors derived from (1) to (8).
Note that, in the deterministic model (10), w is given; θ is the parameter, so it is not decision variable. v D (θ) is the optimal value function parameterized in θ; in other words, it gives the optimal value with any θ in Θ where the parameter set Θ is assumed to be a polyhedron: Set Θ consists of the concerned constraints of decision-makers on θ. Basically, θ includes the power and energy capacities of ESUs, so it should be non-negative; taking the investment cost for example, which relates to the power and energy capacities, the total cost cannot exceed a given budget Π. Nevertheless, Θ can be freely defined.

Ambiguity Set
In reality, the vector w is uncertain due to the inherent natures of renewable power. Combining the scenario-based stochastic optimization and robust optimization, distributionally robust optimization is employed here to optimize the expected renewable energy curtailment, considering the worst probability distribution in a neighborhood of the empirical one built by the representative scenarios.
Given the set W 0 = {w s } S s=1 composed of S representative scenario and their corresponding probabilities {ρ s } S s=1 , we first build the empirical probability distribution, P 0 , by using the following: where δ ρ 0 s is the Dirac measure on ρ 0 s ; thus, P 0 is discrete and supported on W 0 . To construct a neighborhood of P 0 , i.e., the ambiguity set, we need to define the distance between two probability distributions. Given an S-dimensional distribution P supported on W 0 , whose s-th probability is ρ s , the conditions ρ s ∈ [0, 1] and S ∑ s=1 ρ s = 1 are satisfied. The infinite norm-based distance between P and P 0 is as follows [31]: Thus, the ambiguity set based on (13) is developed as follows: where ρ = [ρ 1 , · · · , ρ S ] T and 1 is an all-one column vector with compatible dimensions; ε is an arbitrarily given radius of the ambiguity set, indicating that the infinite norm-based distance between P and P 0 does not exceed ε; meanwhile, a larger ε reflects a higher model conservativeness because some worse distributions are considered. The value of ε can be chosen according to the following [31]: where β is the confidence level; such a choice of ε guarantees the following: Furthermore, M ε I P 0 is a polyhedron and can be expressed as M P 0 = {ρ|Qρ ≤ U}, where Q and U are derived from (14).

Distributionally Robust mp-MILP and Reformulation
If we only consider the empirical distribution, P 0 , the expected renewable energy curtailment, i.e., S ∑ s=1 ρ 0 ĉ Tx s , is minimized. As mentioned above, the distributionally robust formulation minimizes the expectation under the worst distribution in the ambiguity set M ε I P 0 , giving rise to the following DR-mp-MILP: where the vectors with subscript s correspond to the representative scenario w s ; λ is the vector of dual variables. An important feature in (17) is that, given a θ, the feasible set of outer maximization problem is disjoint from that of inner minimization problem; thus, dualizing the outer problem results in the single-level parametric problem shown below: Therefore, to calculate v R (θ) entails solving the mp-MILP (18). For brevity, we aggregate all the coefficients, variables, and constraints in (18), to rewrite it as the canonical mp-MILP:

Solution Methodology
The solution method in this section is based on the mild assumption that the system is feasible given any w s without ESUs (θ = 0), or equivalently, without dispatching ESUs (θ = 0). This implies that (19) is always feasible regardless of y and θ.
Then, we discuss the properties of v R (θ). Binary vector y denotes the charging and discharging pattern of ESUs, and has a finite number of candidates, written as Y = {y i } I i=1 . Given any y i ∈ Y, mp-MILP (19) degenerates into an mp-LP as where v R i (θ) is the corresponding optimal value function; γ is the vector of dual variables. As a result, v R (θ) in mp-MILP (19) can be calculated by comparing v R i (θ), ∀y i ∈ Y in a pointwise minimum way.
Given a θ, consider the dual problem of the mp-LP (20): where Γ = γ A T γ = c, γ ≤ 0 , and b i = b − By i . The duality theory of linear programming indicates that the primal problem (20) and the dual problem (21) share the same optimal value v R i (θ) for any θ; as θ varies in Θ, the optimal solution can always be found at one of the extreme points of Γ. Thus, v R i (θ) can be constructed as follows: Therefore, v R i (θ) is piecewise linear and convex because the point-wise maximum operation preserves convexity; each piece corresponds to a subset of Θ, called critical region (CR). Further, v R (θ) is also piecewise linear but generally not convex.
In Reference [32], a decomposed algorithm is proposed to solve mp-MILP (19), which takes the following four steps:

•
Step 1: Solve MILP by solving MILP (19) by regarding θ as decision variables; the optimal solution is y 1 .

•
Step 3: For the m-th CR with its corresponding piece of v R 1m (θ), build the following MILP, to find a new y that gives a better optimal value, i.e., min x,y,θ c T x where the last constraint excludes the solutions that have been found, and multiplies as the iteration process. The optimal solution of (23), if feasible, is y * 1m , which, along with the m-th CR, solves an mp-LP like (20), whose optimal value function is v R 1m * (θ). Then, the pointwise minimum comparison between v R 1m * (θ) and v R 1m (θ) gives v R 2m (θ), as well as the corresponding CRs.
After checking M 1 CRs, the optimal value function of mp-MILP (19) updates to v R 2 (θ) with M 2 CRs.

•
Step 4: With v R m2 (θ), m = 1, · · · , M 2 and its corresponding m-th CR, go back to Step 3,until (19) finds no new for M K CRs, the algorithm terminates, and the solution of mp-MILP (19) is concluded as follows: where υ k /µ k is the slope/intercept of the k-th piece of v R (θ). Actually, v R (θ) is a pointwise minimum comparison among a group of piecewise linear convex functions, and every single CR in (24) is a convex set, while some CRs may share the same υ k /µ k . Interested readers are recommended Reference [32].

Results
To validate the proposed method, the modified IEEE nine-bus system is utilized, which is sketched in Figure 1. Two generators locate at Bus 1 and Bus 2; one wind farm and one ESU connect to Bus 3. The system data are provided in Table 1. Considering a time horizon of 24 h, i.e., T = 24 h, the load curve is drawn in Figure 2. Ten averagely weighted representative scenarios of wind power are generated by the weather data in Qinghai Province, China, 2019, as in Figure 3. Given β = 95%, the value of ε is set to 0.3, according to (15). Complete data are provided in Reference [34].

Results of Deterministic mp-MILP
Given different scenarios of wind power, (10) turns into a canonical mp-MILP (19), and it can be solved by the algorithm in Section 4. Taking Scenario 3 as an exam the optimal value function is drawn in Figure 4. Furthermore, we define the parameter set as follows:

Results of Deterministic mp-MILP
Given different scenarios of wind power, (10) turns into a canonical mp-MILP like (19), and it can be solved by the algorithm in Section 4. Taking Scenario 3 as an example, the optimal value function is drawn in Figure 4  It can be seen that, without an ESU, the total renewable energy curtailment under Scenario 3 reaches 45.5 MWh. By employing an ESU whose power/energy capacity is higher than 22.75 MW/81.83 MWh, respectively, the curtailment reduces to zero.
The parameter set is divided into three critical regions, named CR1, CR2, and CR3. The polyhedral representations of them are as follows: It can be seen that, without an ESU, the total renewable energy curtailment under Scenario 3 reaches 45.5 MWh. By employing an ESU whose power/energy capacity is higher than 22.75 MW/81.83 MWh, respectively, the curtailment reduces to zero.
The parameter set is divided into three critical regions, named CR1, CR2, and CR3. The polyhedral representations of them are as follows: Then, the corresponding optimal value functions are the following: From (26), we see that the sensitivity of renewable energy curtailment to the capacities of ESU is different as θ changes. When θ falls into CR1, 1 MWh increment of energy capacity will lead to a 0.56 MWh decline of curtailment, while the increase of power capacity makes no sense, implying that the impact of power capacity has saturated; similarly, when θ falls into CR2, a 1 MWh increment of power capacity will cause a 2 MWh decline of curtailment, while raising the energy capacity does not take effect. Moreover, when θ falls into CR3, the curtailment is zero.
Furthermore, the boundary between CR1 and CR2, the red line in Figure 4b, gives valuable information. The points on this boundary are efficient in fully using the capacities of ESUs; from any point on this boundary, increasing power or energy capacity alone will not help reduce the curtailment. Specifically, on this boundary, the ratio of energy capacity to power capacity is 3.60; such a ratio is useful in sizing ESUs in practice.
From an economic perspective, consider the investment cost of ESU as Π = κ P P s + κ E E s , where κ P /κ E is the unit cost of power/energy capacity. Regardless of the value of κ P and κ E , the isocost curves are a cluster of parallel straight lines, as seen from the dashed lines in Figure 5. Combining with Figure 4a, we can find that, given a budget, the parameter on the red line enjoys the minimum renewable energy curtailment, as compared with other points that cost the same.

Results of Distributionally Robust mp-MILP
The optimal value function in mp-MILP (19) is visualized in Figure 6. The parameter set is divided into much more CRs, as compared with the situation in the above deterministic test. As mentioned in Section 4, ( ) is a pointwise minimum comparison among a group of piecewise linear convex functions; every single CR is a convex set, while some CRs share the same / , which are filled with the same color in Figure 6b.
Similarly to Figure 4a, we draw a red line in Figure 6a, on which the points are efficient in fully using the capacities of ESUs from an economic perspective; although this line is not strictly accurate, it does provide a reference for the ratio of energy capacity to power capacity 150 45.49 ⁄ ≈ 3.3 in the sense of distributionally robustness.

Results of Distributionally Robust mp-MILP
The optimal value function in mp-MILP (19) is visualized in Figure 6. The parameter set is divided into much more CRs, as compared with the situation in the above deterministic test. As mentioned in Section 4, v R (θ) is a pointwise minimum comparison among a group of piecewise linear convex functions; every single CR is a convex set, while some CRs share the same υ k /µ k , which are filled with the same color in Figure 6b.
Similarly to Figure 4a, we draw a red line in Figure 6a, on which the points are efficient in fully using the capacities of ESUs from an economic perspective; although this line is not strictly accurate, it does provide a reference for the ratio of energy capacity to power capacity 150/45.49 ≈ 3.3 in the sense of distributionally robustness.
istic test. As mentioned in Section 4, ( ) is a pointwise minimum comparison among a group of piecewise linear convex functions; every single CR is a convex set, while some CRs share the same / , which are filled with the same color in Figure 6b.
Similarly to Figure 4a, we draw a red line in Figure 6a, on which the points are efficient in fully using the capacities of ESUs from an economic perspective; although this line is not strictly accurate, it does provide a reference for the ratio of energy capacity to power capacity 150 45.49 ⁄ ≈ 3.3 in the sense of distributionally robustness.

Results with Different Values of
This part investigates the impact of parameter on the renewable energy curtailment. Parameter reflects the size of ambiguity set and the conservativeness of the distributionally robust mp-MILP. The value of recommended by (15)  According to (15), the value of is mainly determined by , the size of representative scenarios; as the approaches infinity, becomes zero, which means the ambiguity set only includes the empirical probability distribution that is supposed to be exact; meanwhile, the curtailment drops. This concludes the importance of data availability. According to (15), the value of is mainly determined by , the size of representative scenarios; as the approaches infinity, becomes zero, which means the ambiguity set only includes the empirical probability distribution that is supposed to be exact; meanwhile, the curtailment drops. This concludes the importance of data availability.  According to (15), the value of ε is mainly determined by S, the size of representative scenarios; as the S approaches infinity, ε becomes zero, which means the ambiguity set only includes the empirical probability distribution that is supposed to be exact; meanwhile, the curtailment drops. This concludes the importance of data availability.

Results with Different Renewable Energy Penetration Levels
To explore the impact of renewable energy penetration levels on the curtailment, we introduce a ratio r = 0.4, 0.6, 0.8, 1.0, and 1.2 and multiply it with the wind power in every representative scenario; the respective optimal value functions are depicted in

Results with Different Renewable Energy Penetration Levels
To explore the impact of renewable energy penetration levels on the curtailment, we introduce a ratio = 0.4, 0.6, 0.8, 1.0, and 1.2 and multiply it with the wind power in every representative scenario; the respective optimal value functions are depicted in Figure

Discussion
This paper evaluates the impact of energy storage on renewable energy curtailment, using multi-parametric programming; the uncertainty of renewable power is characterized by the ambiguity from following the paradigm of distributionally robust optimization. The original max-min mp-MILP was reformulated into a canonical one by duality theory.
The tests on modified IEEE nine-bus corroborate that the multi-parametric programming framework calculates the explicit relation between ESU capacities and renewable energy curtailment, which enjoys admirable interpretability and convenience for visualization. Specifically, the reference for the ratio of energy capacity to power capacity is concluded; the ratio is about 3~4, with the aim to reduce renewable energy curtailment, alt- More specifically, the optimal value functions can suggest the critical capacities of ESU that can eliminate the curtailment, as θ = [11.1MW, 16.7MWh] for r = 0.4 and θ = [33.3MW, 66.7MWh] for r = 0.6.

Discussion
This paper evaluates the impact of energy storage on renewable energy curtailment, using multi-parametric programming; the uncertainty of renewable power is characterized by the ambiguity from following the paradigm of distributionally robust optimization. The original max-min mp-MILP was reformulated into a canonical one by duality theory.
The tests on modified IEEE nine-bus corroborate that the multi-parametric programming framework calculates the explicit relation between ESU capacities and renewable energy curtailment, which enjoys admirable interpretability and convenience for visualization. Specifically, the reference for the ratio of energy capacity to power capacity is concluded; the ratio is about 3~4, with the aim to reduce renewable energy curtailment, although it is case by case.
Because that the proposed solution method of multi-parametric mixed-integer linear program relies on sampling points from the parameter space, the computational complexity is exponentially related to the size of parameter vector θ; nevertheless, with any fixed θ, it entails solving a linear program, so computational time is generally moderate.
The visualization can be done in, at most, a three-dimensional space. With the growth of the size of θ, it could be difficult in visualization, but we can still derive the parameterized optimal value function, which helps explicitly characterize the impact of parameters.
The results derived from the proposed method, mainly the optimal value function v R (θ), give valuable implications on storage sizing. To find the most effective sizing strategy in reducing the curtailment, the decision-maker can use the following sizing model: To find the most cost-efficient one that means the highest ratio of curtailment reduction to investment cost, use this model: where v R (0) is the curtailment without storage, and γ is the column vector related to the unit investment cost of storage capacity.
To compromise the investment cost and the effectiveness of reducing the curtailment, we resort to the following Nash Bargaining model: Therefore, it is shown that the proposed method is helpful in practice. Furthermore, due to the distributionally robustness, an accurate forecast of renewable power output is not required; only an empirical distribution built by historical data is need. Our model considers all the distributions that are closed to the empirical one, and thus guarantees the performance of evaluation and sizing.