Exact Mixed-Integer Programming Approach for Chance-Constrained Multi-Area Reserve Sizing

An exact algorithm is developed for the chance-constrained multi-area reserve sizing problem in the presence of transmission network constraints. The problem can be cast as a two-stage stochastic mixed integer linear program using sample approximation. Due to the complicated structure of the problem, existing methods attempt to find a feasible solution based on heuristics. Existing mixed-integer algorithms that can be applied directly to a two-stage stochastic program can only address small-scale problems that are not practical. We have found a minimal description of the projection of our problem onto the space of the first-stage variables. This enables us to directly apply more general Integer Programming techniques for mixing sets, that arise in chance-constrained problems. Combining the advantages of the minimal projection and the strengthening reformulation from IP techniques, our method can tackle real-world problems. We specifically consider a case study of the 10-zone Nordic network with 100,000 scenarios where the optimal solution can be found in approximately 5 minutes.

Exact Mixed-Integer Programming Approach for Chance-Constrained Multi-Area Reserve Sizing Jehum Cho and Anthony Papavasiliou , Senior Member, IEEE Abstract-An exact algorithm is developed for the chanceconstrained multi-area reserve sizing problem in the presence of transmission network constraints.The problem can be cast as a two-stage stochastic mixed integer linear program using sample approximation.Due to the complicated structure of the problem, existing methods attempt to find a feasible solution based on heuristics.Existing mixed-integer algorithms that can be applied directly to a two-stage stochastic program can only address smallscale problems that are not practical.We have found a minimal description of the projection of our problem onto the space of the first-stage variables.This enables us to directly apply more general Integer Programming techniques for mixing sets, that arise in chance-constrained problems.Combining the advantages of the minimal projection and the strengthening reformulation from IP techniques, our method can tackle real-world problems.We specifically consider a case study of the 10-zone Nordic network with 100,000 scenarios where the optimal solution can be found in approximately 5 minutes.
Index Terms-Chance constraints, mixed-integer programming, multi-area reserve sizing, probabilistic constraints.

I. INTRODUCTION
I N EUROPE, transmission system operators (TSOs) are in- creasingly coordinating their system operations in response to the pan-European coupling of electricity markets [1].One of the objectives of this coupling is to organize a system that encompasses multiple areas for dispatching balancing energy from frequency restoration reserves in real time or close to real time (the MARI and PICASSO platforms). 1 An important problem of interest that is emerging as a result of cross-zonal coordination in balancing is to allocate the right quantities of reserves in the right locations of the network while accounting for possible congestion in the transmission network.This problem is referred to as reserve sizing or reserve dimensioning, with the associated challenge of reserve deliverability [2], [3], depending on the context.
Article 157 of the System Operation Guideline (SOGL) of the European Union [4] explicitly specifies probabilistic requirements for reserve sizing.The Nordic System Operation Agreement (SOA) [5] is an example of an effort for the coordinated operation of frequency reserves among the Nordic countries in response to the SOGL.A recent ENTSO-E report [6] by Danish TSO Energinet demonstrates the continual pursuit in the direction of multi-area reserve sizing in accordance with article 157 of the SOGL.
There exist a number of papers that attempt to address the multi-area reserve sizing problem in the literature.However, much of this literature [7], [8], [9], [10], [11] focuses on either chance constraints or transmission constraints without treating these aspects jointly.Other literature [12], [13] considers these two aspects simultaneously; nonetheless, the underlying probabilistic distributions are assumed to belong to a specific class.Recent literature [14], [15], which has been proposed by the authors of this article, accounts for both of these characteristics.In [14], the authors define a chance-constrained formulation for the problem, but suggest a heuristic method that is not guaranteed to furnish the optimal solution.Subsequently, in [15], the authors attempt to solve this two-stage mixed-integer programming to optimality by applying integer programming techniques.However, the method proposed by the authors is not scalable to the size of realistic problems.
Starting from [16], there exists a strand of literature about using integer programming techniques for solving chanceconstrained optimization problems [17], [18], [19], [20], [21], [22], [23], [24], [25], [26].However, not all of these methods can be directly applied to our problem due to the strict requirements for the structure of the problem that they impose.One of the common characteristics of these techniques is that they all try to reduce the gap between the linear programming relaxation solution and the true integer optimal solution in order to enable Branch-and-Bound algorithms to solve the problems more efficiently.These 'gap tightening' techniques are typically designed for a specific structure of the problems.Unfortunately, minor changes to the underlying models can render these techniques no longer applicable.
Our problem is formulated as a two-stage chance-constrained problem, because we decide reserves in the first stage and react in the second stage given the reserves that we have decided.Only some of the techniques referenced previously [19], [20], [22] can be directly applied to this form of the problem.However, these techniques do not scale well in large-scale instances of two-stage problems.For example, [15] is based on the method proposed in [19]; however this method is not scalable.This is because the methods need to scan all scenarios and solve a linear program for each one every time the algorithm generates a single inequality to tighten the gap.For large-size instances, the number of scenarios and the number of generated inequalities are too high.
Nevertheless, by reformulating the two-stage model to an one-stage problem, other methods can be utilized, such as [16], [17], [18], [24], [25], [26].These methods are designed for one-stage joint2 chance-constrained problems.Thanks to the simpler structure that they target, these methods are more powerful and can deal with larger instances.However, the reformulation from a two-stage problem to an one-stage problem is not straightforward.In general, finding a projection from a space with a higher dimension is difficult.Applying general-purpose projection methods [27], [28] requires a prohibitive amount of computation, especially for large-scale instances.Additionally, finding a minimal representation for this projection is important because performance in joint chance-constrained programs depends significantly on the number of constraints that are contained in the probabilistic constraints.In this article, we provide a closed-form minimal projection formulation.This allows us to develop a method that is scalable to realistic problems, as we demonstrate in a case study at the end of the article.
Notice that, although it is a necessary step for the whole approach, having a minimal projection formulation itself is not enough to solve instances of realistic scale to optimality.Tightening through integer programming techniques is also an essential part.The projection enables us to utilize more powerful techniques, but it does not solve the problems in itself.Among several integer programming methods that are available for joint chance-constrained programs, we introduce a method that is powerful but also simple to implement [16].
In Section II, we present a formal definition of the multi-area reserve sizing problem.We introduce our minimal projection formulation with mathematical proofs.In Section III, we further reformulate it in a form where we can use a commercial solver and introduce integer programming techniques for tightening the LP relaxation gap, which is crucial for the performance of the method.First, we start with the sample approximation scheme to tackle the probabilistic constraint.Then, we present integer programming techniques based on mixing inequalities and a strong extended formulation using this concept.In Section IV, we formally present our strengthened minimal projection formulation, which is a combination of the minimal projection formulation in Section II and the integer programming techniques introduced in Section III.Detailed explanations for implementing our algorithm are provided in this section.In Section V, a realistic case study of the Nordic system is presented.We compare our method with the heuristic from [14] and another alternative exact method from [15].In-depth simulation results over different sample sizes are also presented.Section VI summarizes the analysis and proposes areas of future research.

II. PROBLEM FORMULATION
For a network G(Z, E), let r + z [resp.r − z ] denote the size of upward [resp.downward] balancing capacity of reserve for each zone z.Our goal is to minimize the sum of r + z and r − z for all the zones in G(Z, E).Note that the objective function can be extended straightforwardly to the case where total procurement costs are considered through balancing capacity offers.In this case, the coefficients of r +/− z would be different values from +1, but the method in this article can manage this type of adjustment.F + [resp.F − ] denotes the feasible region for r + [resp.r − ] representing the region where the capacity of reserve can cover imbalances δ z for each zone z in the network G(Z, E).Formally, F +/− are defined as (1) and (2).
Here, p z and f e are the amounts of balancing energy activated at zone z and the flow from z 1 to z 2 , where e = (z 1 , z 2 ), given that link e has capacity limits T + e and T − e in the reference and opposite direction respectively.The equations in the first lines denote the power balance equations for each zone z.The inequalities in the second line impose that the activation of balancing energy cannot exceed the amount of available reserve.Flow limits are imposed in the last inequalities.Note that the values of δ z , T +/− e can vary under different scenarios, as we discuss in the sequel.
Notice that, in this article, we assume that power flow constraints are approximated using a transportation network model.This assumption is aligned with the fact that, within MARI, the platform for the activation of manual frequency restoration reserve, the network will be approximated using an ATC (Available Transfer Capacity) transportation-based model [29], [30] at the launch of the platform.
This is a two-stage chance-constrained formulation where the first-stage variables are r + , r − and the second-stage variables are p z and f e .There are two different directions that we consider for reformulating the constraints (3b).In this article, we use a projection method to represent the feasible regions of the firststage variables r + , r − explicitly in the space of the first-stage variables.Notice that it is more common to use the second-stage variables p z and f e in order to represent r +/− ∈ F +/− .One way to use this approach with a heuristic method will be introduced and compared with our method in Section V-A.

A. Minimal Projection Formulation
In this subsection, we introduce a way to reformulate the constraints (3b) so that they represent r +/− ∈ F +/− explicitly in the space of r +/− .Intuitively, one might surmise that imposing that the reserves should be greater than or equal to the lowest possible aggregate imbalances corresponding to each element of the power set 3 of the vertex set would suffice to represent F +/− explicitly with r +/− .Although this is not trivial, it happens to be correct.However, this representation is not minimal, since there are inequalities that can be represented by a linear combination of other inequalities.In this section, we present a compact representation that is minimal.It turns out that using only the subsets of the vertex set whose elements are "connected" on the network is enough to obtain a minimal projection.We, therefore, define this set as a connected vertex set in order to develop our analysis.For each element of the connected vertex set, imposing that the sum of the reserve capacities of the zones included in the element is greater than equal to the size of the total imbalances in the element minus the maximum input (or output) flow of the element, is sufficient to represent F +/− .This statement can be formally stated as Theorem 2.1, and we prove it in the sequel.We commence by formally defining the following objects: connected vertex set, and maximum input/output flow.
Definition 2.1 (Connected Vertex Set): For a graph G(V, E), the connected vertex set W(G) is defined as follows: (4) where V (P ) denotes the set of vertices in the path P .
Example 1 (Connected Vertex Set): For the graph in Fig. 1 , {1, 2, 3} is an element of a connected vertex set, whereas {1, 4} is not.For a vertex set {1, 2, 3}, as an example, when v = 2, w = 3 there exists a path . This is true for all possible combinations of v, w ∈ {1, 2, 3}, so this is an element of the connected vertex set.On the other hand, for a vertex set {1, 4}, when v = 1, w = 4, there is no such path whose vertex sets are subsets of {1, 4}, since all the paths between v = 1 and w = 4 contain either 2 or 3, which are not elements of {1, 4}.
Intuitively, if all the elements in a vertex set have an edge that connects these vertices to any of the other elements in the vertex set, such vertices are elements of a connected vertex set; hence 3 Recall that the power set of a set is the set of all subsets of the set.the name of the definition.For the 5-zone example in Fig. 1, the connected vertex set is Notice that maximum input/output flows are properties of the network since they are defined only on the basis of line capacities Example 2 (Maximum Input/Output Flow): In Fig. 2 , every edge of the directed graph has flow capacities in both directions.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
For the edge (1,2), the maximum capacity for the direction , and the maximum capacity for the direction 2 → 1 is T − (1,2) .Maximum input/output flows are defined on an edge subset E = {(1, 2), (2,4), (4,3)} ⊂ E for a vertex subset S. In Fig. 2, S = {2, 4}, so the maximum input flow is 3) , and the maximum output flow is ) .Let us return to our problem.For the sake of brevity, let us denote W(G(Z, E)) for the network G(Z, E)) for our original problem as W(G) from now on.Using this new notation, we define the three following sets F, F p , F r .
Notice that F is a set defined with the same types of constraints as F +/− from ( 1) and ( 2), but considering upper and lower bounds for the balancing energy p z at the same time in (8).Another difference is that F is in the space of (r +/− , p, f).The goal is to find a projection onto the space of r +/− only.The resulting projection is the set F r and F r = F + ∩ F − .F p is a projection of F onto the space (r +/− , p), and it is used as an intermediate step to go from F to F r in order to prove that F r is indeed a projection of F .Since the proof for the set F is more general than the proofs for (1) or ( 2), we present the proof for F in this section.From now on, for a set A defined in the space of variables (x, y), we denote P roj (x) (A) as the projection of the set A onto the space of x.Theorem 2.1: The proof is based on two steps.First, in Claim 2.1.2,we show that the projection of F onto the space of (r + , r − , p) is F p .Second, in Claim 2.1.1,we show that the projection of F p onto the space of (r + , r − ) is F r .Claim 2.1.2and Claim 2.1.1 together imply that P roj (r The proofs for Claims 2.1.1 and 2.1.2are provided in the appendices.
Corollary 2.2: , where Theorem 2.1 and Corollary 2.2 show that the set F r is indeed an explicit representation of the projection of F on the space of the first-stage variables r +/− , resulting in F +/− = F +/− r .Theorem 2.3: F r is a minimal representation on the space of (r + , r − ).
Proof: Since the proof for the set of inequalities ( 12) is similar to the case for (13), we show here the case for (12).In order to show that ( 12) is a minimal representation on the space of r − , to arrive at a contradiction, first let us assume that there exists a set S ∈ W(G) such that there exist mutually different sets S 1 , . . ., S n ∈ W(G) by which the inequality constructed in the form of (12) dominates the inequality for the set S .Formally, this means that there exist coefficients α 1 , . . ., α n ≥ 0 that satisfy the following conditions ( 14) and ( 15): In order to satisfy the inequalities ( 14) and ( 15) for all possible values of r − v and δ v , Σ i:v∈S i α i = 1 for all v ∈ S and Σ i:v∈S i α i = 0 for all v ∈ V \ S .This implies that for all i ∈ {1, . . ., n}, S i ⊆ S and n i=1 S i = S .Notice that the left-hand side and the right-hand side of ( 14) are equal, and (15) becomes Since the right-hand side of ( 16) where this contradicts the initial assumption.Q.E.D. Corollary 2.4: F +/− r is a minimal representation on the space of r +/− .Theorem 2.3 and Corollary 2.4 show that the sets of inequalities ( 12) and ( 13) are indeed minimal representations for the projection of F (or F +/− ).Notice that the inequalities (12) and ( 13) have an intuitive explanation.For certain combinations of zones, the sum of the reserve capacities should cover the total imbalances for the zones net of the maximum input [resp.output] flow for upward [resp.downward] reserves.For the case of infinite line capacities, where T +/− is infinity, our multi-area problem amounts to a single-zone problem.This can be checked from ( 12) and ( 13), where the only non-redundant constraints are when S ∈ W(G) is equal to the set of all the zones Z, which Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. is equivalent to the case where we aggregate all the zones in one region.
Example 3 (Minimal Projection): In order to illustrate the reformulation more explicitly, we provide a three-node example in Fig. 3. Firstly, the connected vertex set for the three-node graph is Accordingly, we write down all the constraints for F + r as follows: r can be written down in a similar fashion.
As in Example 3, the representation r +/− ∈ F +/− in constraint (3b) can be replaced by A +/− r +/− ≥ ξ +/− , where A + and A − are linear maps, and there are uncertain parameters only in the right-hand-sides ξ + and ξ − .The resulting new formulation is only involving the first-stage variables.Now that we have derived an explicit representation with inequalities, this new formulation enables us to apply the integer programming techniques for joint chance-constrained programs, which we address in the next section.Thanks to Corollary 2.4, we are guaranteed that our representation is minimal, so that the number of inequalities within the probabilistic constraints is minimized.This is crucial for the computational performance of our method, since, as the number of inequalities increases, the linear programming relaxation gap tends to be larger.We can now introduce a tightening step using integer programming techniques.In the following section, we introduce one such technique known as strong extended formulation.

III. STRENGTHENED FORMULATION
Thanks to the projection formulation in the previous section, our problem (3) can be represented as follows. min The main difference between this representation and (3b) is that the polyhedron on the space of the variables r +/− in the probabilistic constraint is now explicitly known in the form of a set of inequalities.

A. Sample Approximation
Sample approximation is a standard way of dealing with probabilistic constraints in chance-constrained problems.This approach is first introduced in [31].As the name suggests, it approximates a probabilistic constraint by sampling several scenarios for uncertain parameters.By doing so, we do not need to restrict the underlying uncertainty distribution to be of a specific type, such as Gaussian.Assuming that we have a sufficiently large sample set, this method is guaranteed to converge to the true optimal solution of the optimization problem.Since it is often the case that we do not know the underlying distribution for the uncertain parameters, but do have access to historical data, this approach is widely used in many applications. Given where u +/− i = 1 indicates that, under scenario i, the constraint of balancing the system is violated; thus, the number of violated scenarios should be less than or equal to +/− N , according to the direction of the reserves.
The reformulation of ( 18) is a mixed-integer linear programming problem.Although it is a form that can be plugged into a commercial solver, this step alone cannot solve large-scale instances to optimality due to the big LP relaxation gap.In order to close this gap, we need to exploit the specific structure of the formulation.In (18), our main constraint has a form that is widely studied, the so-called mixing set.This allows us access to a set of techniques that help us close the LP relaxation gap.In the following sub-sections, we briefly introduce the mixing set and how we utilize this set for closing the LP relaxation gap, followed by the technique we use for our method.

B. Mixing Set and Mixing Inequalities
Given N scalars h i for i ∈ [N ], a mixing set is defined as (19) Note that, when y is equal to the j-th row of A +/− r +/− , and h i = ξ ij , where ξ ij denotes the j-th component of the vector ξ i , our formulation in (18) has the form of a mixing set.In fact, there are multiple mixing sets, as many as the number of rows Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of A +/− r +/− .Each of these sets is tightened independently through the special types of inequalities named mixing inequalities.
Mixing inequalities, also known as star inequalities, are developed by Atamturk et al. [32] and Gunluk and Pochet [33].If we assume that h 1 ≥ h 2 ≥ • • • ≥ h N without loss of generality, the mixing inequalities for the mixing set (19) are defined as where t 1 < • • • < t l and h t l+1 := 0. It is known that the mixing inequalities are valid and are sufficient for defining the convex hull of the mixing set (19).Defining a convex hull with a certain set of inequalities implies that the LP relaxation gap is zero; in other words, we can solve the Integer Programming problem by solving a Linear Programming problem with these convex hull defining inequalities.This fact can be used even for more complex problems in which we are not necessarily able to characterize the convex hull of the problem.When problems are more complex, we can no longer guarantee a zero LP relaxation gap; however, these mixing inequalities are still valid for those problems, and this aids significantly in reducing the LP relaxation gap, resulting in better performance for solving the original Integer Programming problem with binary variables u.
Furthermore, for the following set with a cardinality constraint induced by a reliability criterion , where we define q = N , the strengthened mixing inequalities with t 1 < • • • < t l and h t l+1 := h t q+1 are facet-defining for conv(G) if and only if t 1 = 1 [16].An inequality is facetdefining for conv(G) when the inequality is crucial for defining the convex hull of the set G. It implies that these inequalities are expected to be the most effective for reducing the LP relaxation gap among similar types of inequalities, since the convex hull is the smallest convex set containing all the feasible solutions.
In this article, we assume that all the scenarios have equal probabilities.There are similar techniques for knapsack constraints for unequal probabilities.In this case, the inequalities are less tight.However, there is no significant loss of generality since we can approximate the problem with equal probabilities by resampling from the original distribution.The readers who are interested are referred to [16].

C. A Strong Extended Formulation
Although these inequalities are useful for tightening the LP relaxation gap, there are additional steps required to use them.One might first consider completely enumerating all such inequalities.However, since there are exponentially many such inequalities, this will not be a practical approach.Alternatively, there are ways to add a subset of the inequalities on-the-fly, while we are running the Branch-and-Bound algorithm.This is a potential way forward; however, it is not straightforward to implement such algorithms.In this article, we introduce a method that is very easy to implement and achieves excellent performance for our application.
The essence of this method is to use a very compact formulation that has the same effect as when we add all the inequalities (22) by introducing a new set of variables.This new formulation is called a strong extended formulation and is first introduced in [16].The name 'extended' comes from the fact that there are additional variables to the original formulation, and 'strong' refers to the fact that this new formulation is as strong as adding the entire exponential family of valid inequalities (22) to the original one.
Formally, the extended formulation of G in ( 21) is defined as follows: where Theorem 3.1: (Theorem 6 from [16]) P roj (y,u) (EG) = G.Moreover, the projection of the linear relaxation of EG is the linear relaxation of G with all the inequalities (22) added.
As Theorem 3.1 shows us, by simply introducing q binary variables w, we can tighten the LP relaxation gap.Notice that (24) are easily implemented in a commercial solver.In the next section, we present the end result of the combination of Sections II and III.The minimal projection formulation in Section II is re-formulated with the idea of the strong extended formulation in Section III.

A. Formulation
The strengthened minimal projection formulation using F r (12), (13) and EG in ( 23) is formally defined as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
For S ∈ W(G), where σ +/− S,i are the permutations that rearrange the indices as

B. Implementation
The formulation of ( 25) is a Mixed-Integer Linear Programming problem which can be directly solved by commercial solvers such as CPLEX and GUROBI.However, in order to implement the algorithm, three elements are required: the connected vertex set W(G), the coefficients of mixing sets h The W(G) Generation algorithm presented in this article has a worst-case complexity 4

of O(|E| • 2 |V |
).This is significantly lower than the complexity of general-purpose projection methods.For instance, the Fourier-Motzkin-Elimination algorithm [27], which is a well-known general projection method, is known for its double exponential complexity, that is O((|E| + |V |) 2 (|E|+|V |) ) for our case.A recently developed general algorithm [28] enjoys a single exponential complexity, but with much higher exponents; namely O((|E| 3 ) for our problem.
2) Sorting: Now that we have |W(G)|, we are ready to generate the remaining elements h +/− S,i and σ +/− S,i .An easy way to obtain the parameters h S,i is to calculate the righthand-side for each scenario i in (26) and sort these right-hand side parameters in non-increasing order.Then, the resulting non-increasing sequence becomes h +/− S,i and the corresponding permutation of the indices becomes σ +/− S,i .As the size of the connected vertex set can be exponential with respect to the number of zones, the time for the process of calculating h +/− S,i and σ +/− S,i is non-negligible.In the next section, we present the calculation time for pre-processing as well as the solver time for the optimization problem when we present the computational results for a case study.

V. COMPUTATIONAL RESULTS
In this section, we compare our strengthened minimal projection method with several alternatives.First, we introduce another way to reformulate the problem using the so-called "Big-M" method.Even though this approach is not scalable in practice, it is also a basis of a heuristic method introduced in [14].We compare both ways of using this reformulation with our method.Additionally, we also compare with another exact method using Benders' Decomposition.

A. Big-M Based Formulation With Second-Stage Variables
In this alternative formulation, we also use the sample approximation approach to deal with the probabilistic constraints of (3b).
Similarly to (18), binary variables u +/− are introduced to represent whether the probabilistic constraint is violated or not under scenario i.The difference here is to use the second-stage variables p and f directly and use a Big-M method to represent the logical expression 5 u . This is achieved by introducing slack variables l +/− .Slack variables l +/− zi are non-zero only when u +/− i = 1, enabling the power balance constraints to be violated.Notice that l +/− zi is bounded by max{0, −δ zi } or max{0, δ zi }; they are the upper and lower bounds for p zi in absolute value.Thanks to this large bound for l +/− zi , when u +/− i = 1, there exists a feasible solution with p i = 0, f i = 0, which allows scenario i to not be accounted for when calculating the size of reserves r +/− .Unfortunately, it is well known that formulations using the Big-M method are not practical for solving problems to optimality due to large LP relaxation gaps.However, the model of (29) can be used for developing a heuristic method in order to find a feasible solution.In [14], for example, the authors first solve the LP relaxation of (29) and fix u +/− i to 1 for the indices in which the optimal solutions for the LP relaxation u * +/− i are in the sets of the +/− N largest values when the solutions u * +/− i are sorted in descending order.

B. Comparison With an Exact Method
In this subsection, we compare our method with another exact method that also guarantees to solve the problem to optimality [15].This article uses the theory in [19], which is based on Bender's Decomposition [34].Roughly, this method tries to solve the two-stage chance constraint problem directly instead of reformulating it.It does so by generating inequalities of the type of ( 22) through Bender's Decomposition.The computational results of this alternative method are much faster than solving the Big-M formulation of (29) to optimality.However, as discussed in the introduction, this direct approach to solving twostage chance-constrained problems is not scalable.The method 5 We denote F +/− i as the sets F +/− where all the random variables (δ z , ) are replaced by their realizations for scenario i; namely (δ zi , T in [15], for example, requires approximately 30 minutes to solve instances of four zones with 5,000 samples to optimality.On the other hand, our method using the formulation of (25) can achieve the optimal solution for the same size of instances within 1 s.One of the reasons for the subpar performance of this alternative method is the inequality generation step.Methods based on Bender's Decomposition in chance-constrained problems, in order to generate a single inequality that is similar to (22), a linear program should be solved N times, the number of scenarios.As the size of the problem increases, the number of inequalities and scenarios increase at the same time.On the contrary, through the minimal projection step and the additional strengthening step, our formulation is already in a compact form that does not require any inequality generating steps.Since this method has been found to not scale to large instances, we proceed to compare our method with a heuristic method using the Big-M formulation of (29).

C. Case Study: Comparison With a Heuristic Method
For the comparison, a case study of the Nordic system is considered.In this case study, as indicated in Fig. 4, three Nordic countries (Norway, Sweden and Finland) are involved, and they account for 10 bidding zones with 15 links.The reference data for imbalances for each zone and the network capacity are sourced from [35].For the imbalances, we generate samples from a normal distribution with zero mean and a standard deviation equal to the reference imbalances.For the network capacity, we add perturbations to the reference data for each sample.The perturbations are distributed according to a normal 6.Sensitivity analysis for the sample approximation approach over sample size when epsilon is 1%.distribution with zero mean and a standard deviation equal to 5% of the value of the reference data .For all the Figs.5-7, the bar charts refer to the mean of 100 simulations in which the middle lines indicate the standard deviation.Throughout the case study, we use GUROBI version 9.51 as optimization solver with JuMP embedded in the programming language Julia, and the computing equipment with a SkyLake CPU (2.3 GHz).
We compare our method in Section IV and the LP based heuristic method in [14], introduced in Section V-A.In Fig. 5, the results comparing (a) the optimal reserve sizes and (b) the total solving time6 are presented for varying degrees of reliability levels ( ) when the sample size is N = 25, 000.The Minimal Projection Method can be solved notably faster than the LP Based Heuristic, and finds the optimal solution.This seemingly counter-intuitive result can be explained by the fact that the strengthened minimal projection formulation (25) often has a smaller size than the formulation (29) in terms of the number of variables and constraints.Notice that q << N, thus each set of constraints in (25) is repeated q times whereas that in (29) is N times.Additionally, even though the size of the connected vertex set W(G) is exponential, it is often the case that h +/− in (26) are all negative, resulting in adding redundant constraints that can be ignored or automatically removed in pre-processing steps of commercial optimization solvers.This phenomenon happens more often when the capacities of lines T +/− are sufficiently large compared to the level of imbalances δ.In an extreme case where T +/− is infinity, the only constraints left are z∈Z r +/− ≥ ∓ z∈Z δ z , which is equivalent to a single-zone problem.
The gap between the optimal solution and the sub-optimal solution from the LP method when = 1% is around 18.9%.As becomes smaller, since q also becomes smaller, the optimization problem becomes less complex, resulting in a smaller gap between the optimal solution and the sub-optimal solution and a faster solving time.However, the high value of the standard deviation in the optimal solution for = 0.1% implies that the sample size is not sufficient.When we increase the sample size, then the gap also increases.
Additionally, in Fig. 6, we present a sensitivity analysis of the sample approximation approach with respect to the sample size for our case study.There are two issues when the number of samples is low.First, the variance of the optimal objective function value is high.In Fig. 6, this is captured with the coefficient of variance (CV), which is the variance divided by its mean.One can observe that CV is the highest (1.1%) for the case of N = 10, 000, and as the sample size increases, the CV decreases and the value of reserve stabilizes.Secondly, a typical phenomenon of the sample approximation method that is introduced in [31] is an underestimation of the true objective function value, in the sense that the resulting optimal objective function value tends to be lower than the true optimal objective function value of the problem.This can be observed in Fig. 6.
Lastly, we analyze the solving time for the Strengthened Minimal Projection Method over different sample sizes in Fig. 7. Notice that = 1% is the most computationally complex problem among the three different levels of .Our method allows us to solve N = 500, 000 samples in less than 30 minutes in terms of optimization time.In general, the data pre-processing time is non-negligible due to the exponential size of the connected vertex set W(G); however, the bottleneck complexity is O(N logN ) due to the sorting algorithm that is still scalable.In practice, if one needs to solve the problem dynamically, adding new samples to an already-sorted list (which can be expected to be the case in practice, based on information communicated to us by Nordic TSOs) is much easier than sorting the entire list, and in this case the data pre-processing time is negligible.

VI. CONCLUSION
In this article, we propose a novel method for solving the chance-constrained multi-area reserve sizing problem to optimality.After identifying a minimal representation of the projected set of our feasible region, we use integer programming methods to strengthen our formulation.This approach can deal with instances of realistic size, and this is shown in a case study of the Nordic system.Since a transportation-based network is used for our method, this approach can also be used in different domains.
In future work, it is possible to extend the model with different approximations of power flow constraints.A DC (Direct Current) approximation can be an option.From the perspective of better calculation, we can apply more recent IP techniques for solving the minimal projection formulation.

Algorithm 2: Finding
the case where R 1 = ∅.The lower bound of (A.8) ≤ the upper bound of (A.7) is implied by (13) and the upper bound of (A.9) ≥ the lower bound of (A.7) is implied by (12).For showing why the lower bound of (A.8) ≤ the upper bound of (A.9), pick 12) for S 2 \ S 1 and (13) for S 1 \ S 2 7 using Lemma A.1, which is equivalent to the lower bound of (A.8) for S 1 ≤ the upper bound of (A.9) for S 2 .Thus, pv 1 satisfying (A.7)-(A.9)exists for the case where R 1 = ∅.
For the next step of mathematical induction, assume that, for i ≥ 1, there exists pv k for 1 ≤ k ≤ i satisfying (A.7)-(A.9).For R i+1 = R i ∪ v i and v i+1 ∈ V \ R i+1 , our goal is to show that all the possible combinations of the upper bounds and the lower bounds from (A.7)-(A.9)can be implied by other inequalities so that we can show that pv i+1 exists.First, we show it for the combinations of upper bounds and lower bounds between (A.7) and (A.8)-(A.9).Here, we show one out of the two cases: the lower bound of (A.8) ≤ the upper bound of (A.7).The other case can be 7 It is possible that shown in a similar fashion.The set W(G) can be divided into two cases : i) R i+1 ∩ S = ∅ and ii) R i+1 ∩ S = ∅.For the case i), Σ w∈R i+1 ∩S pw = 0 and Σ w∈S\{R i+1 ∪v i+1 } r+ w = Σ w∈S\v i+1 r+ w , so (13) implies the lower bound of (A.8) ≤ the upper bound of (A.7).For the case ii), from the set {v : v ∈ R i+1 ∩ S}, pick the node with the largest index l.Observe that Σ w∈R i+1 ∩S pw = Σ w∈R l ∩S pw + pv l and Σ w∈S\R i+1 r+ w = Σ w∈S\{R l ∪v l } r+ w .This can be proven by contradiction.Assume that it is not true.Then ∃v m such that m = l, v m ∈ R i+1 , v m ∈ R l , and v m ∈ S.This contradicts the fact that l is the largest index.Thus, (A.8) with R l and v l implies the lower bound of (A.8) ≤ the upper bound of (A.7).
For showing why the lower bound of (A.8) ≤ the upper bound of (A.9), pick Since it is similar in the other cases, here we only show the argument for the case ii) where  12) for S 2 \ S 1 using Lemma A.1, following a similar process as in (A.11) and (A.12) we get the inequality, which is equivalent to the lower bound of (A.8) for S 1 ≤ the upper bound of (A.9) for S 2 .Thus, pv i+1 satisfying (A.Proof: Since it follows an almost identical reasoning, we only show the case of Maximum Output Flow.For the sake of compactness, without loss of generality, we leave out the conditions (v, w) ∈ E or (w, v) ∈ E under the summation sign.Notice that O(S|E ) consists of the terms related to T + (v,w) and those of T − (v,w) .In this proof, the patterns for T + (v,w) and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
T − (v,w) are exactly same and what is important is the relationship of summations, so we omit T + (v,w) and T − (v,w) over the course of the equations.Notice that the right-hand-side can be written as follows: O(S  Proof: For the sake of compactness, without loss of generality, we leave out the conditions (v, w) ∈ E under the summation sign.Proof: First, we show that P roj (r + ,r − ,p) (F ) ⊆ F p .Notice that (8) and (11) are identical.So, it suffices to show that (7) and ( 9) imply (10).From (7),  9) imply (10).Second, we show that F p ⊆ P roj (r + ,r − ,p) (F ).It suffices to show that for all (r + , r− , p) ∈ F p , there exists f such that (r + , r− , p, f ) ∈ F. We

Fig. 1 .
Fig. 1.A graph with 5 zones for illustrating the definition of a connected vertex set.
and the size of the connected vertex set |W(G)| is 21.Definition 2.2 (Maximum Input/Output Flow): For a directed graph G(V, E) where ∀e ∈ E, f (e) denotes the flow in e and −T − e ≤ f (e) ≤ T + e , for all S ⊆ V, E ⊆ E, the maximum input flow I(S|E ) and the maximum output flow O(S|E ) on E are defined as follows:
+/− S,i and the permutations σ +/− S,i for all S ∈ W(G) and i ∈ [N ]. 1) Generation of Connected Vertex Set: The connected vertex set of a graph G = (V, E) can be generated using Algorithm 1.The size of the resulting connected vertex set W(G) varies according to the topology of the graph G.When the graph is radial and all the nodes are connected (such as chains), the size of the connected vertex set is minimal and it is |W(G)| = 1/2 • |V |(|V | + 1).The worst-case scenario is when the graph is a complete graph where all the nodes are connected to each other.In this case, the size of the connected vertex set is |W(G)| = 2 |V | − 1.For the graph in Fig. 1, there are 5 zones, so the size (|W(G)| = 21) is within the range of 4 • 5/2 = 10 ≤ 21 ≤ 2 5 − 1 = 31.

Fig. 4 .
Fig. 4. Bidding zones and transmission network lines for a case study of the Nordic countries.

Fig. 5 .
Fig. 5. Comparison between Strengthened Minimal Projection Method and the LP Based Heuristic Method when the sample size is N = 25, 000 in terms of (a) optimal objective function value and (b) total solving time.

Fig. 7 .
Fig. 7. Solving time for the Strengthened Minimal Projection Method over sample size when epsilon is 1%.
and we can get the same results as (A.11) by summing up(12) or(13) for S A and the same for S B .