Pool & Discard Algorithm for Chance Constrained Optimization Problems

In this paper, we describe an effective algorithm for handling chance constrained optimization problems, called the Pool & Discard algorithm. The algorithm utilizes the scenario approximation framework for chance constrained optimization problems, and the warm-start and problem modification features of modern solvers. The exploitation of the problem structure and efficient implementation allows us to considerably speed up the computations, especially for large instances, when compared with conventional methods.


I. INTRODUCTION
This article describes a novel method for handling chance constrained optimization problems that was developed in the author's dissertation [1]. The introduction into the topic of chance constrained optimization is derived (more or less directly) from [2] -with most of the used notation adapted from [2] as well. Let X ⊆ n x be a convex and closed domain of optimization and consider a family of constraints x ∈ X ξ parameterized in ξ ∈ . The uncertain parameter ξ describes different instances of an uncertain optimization scenario. We adopt a probabilistic description of uncertainty and suppose that the support for ξ is endowed with a σ -algebra D and that a probability measure P is defined over D. The probability measure P describes the probability with which the uncertain parameter ξ takes value in . Then, a chance constrained optimization program is written as: Here, we assume that the σ -algebra D is large enough, so that {ξ : x ∈ X ξ } ∈ D, i.e. {ξ : x ∈ X ξ } is a measurable set. Also, linearity of the objective function can be assumed without The associate editor coordinating the review of this manuscript and approving it for publication was Sun-Yuan Hsieh . where y is a scalar variable.
In the CCP (1), constraint violation is tolerated, but the violated constraint set must be no larger than . The parameter allows us to trade robustness (in terms of the probability of constraint violation) for performance (in terms of the optimal objective value): the optimal objective value J * of CCP is a decreasing function of and provides a quantification of such a trade-off. Depending on the particular application (the range of applications is quite wide), can take different values and has not necessarily to be thought of as a ''small'' parameter.
Chance constrained programming has been around for more than half a century, at least since the work of Charnes, Cooper and Symonds in the fifties, see [3]. In [3], however, only individual chance constraints were considered. Joint probabilistic constraints, as in (1), were first considered by Miller and Wagner, [4], in an independent context, while a general theory is due to Prékopa, see [5], [6]. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Prékopa was also the one to introduce the convexity theory based on logconcavity, which was a fundamental step toward solvability of a large class of chance constrained problems. The books [7] and [8] provide an excellent and broad overview on logconcavity theory in stochastic programming, and related results. Yet another study about the convexity of chance constrained problems is [9], while convex approximations of chance constrained problems are considered in [10], [11], and [12]. Stability of the solution under perturbation of the chance constrained problem is studied in [13] and [14]. Although chance constrained problems can be efficiently solved in some special cases, it remains true that the feasible set of CCP is in general non-convex in spite of the convexity of the sets X ξ . Therefore, an exact numerical solution of CCP is, at least in general, extremely hard to find.

II. SAMPLE COUNTERPART APPROACH
We can view the variable x ∈ X ⊆ n x as the ''design variable''. The family of possible instances is parameterized by an ''uncertainty vector'' ξ ∈ ⊆ n ξ . Then, the prototype optimization problem consists in minimizing a linear objective c T x, subject to that x satisfies the constraints is a scalar-valued function that specifies the constraints. Note that considering scalar-valued constraint functions can be assumed without loss of generality, since multiple constraints g 1 (x, ξ ) ≤ 0, . . . , g m (x, ξ ) ≤ 0 can be expressed by a single scalar-valued constraint by the position g(x, ξ ) = max i=1,...,m g i (x, ξ ). Although convexity is preserved by this operation, other valuable properties, such as linearity or differentiability, are lost. In typical situations, has infinite cardinality, i.e., it contains an infinite number of possible instances for ξ . Assumption 2.1 (Convexity): For each ξ ∈ the sets X ξ are convex and closed. Assumption 2.1 requires convexity only with respect to the design variable x, while generic nonlinear dependence with respect to ξ is allowed.
Depending on the situation at hand, the measure P can have different interpretations. On one hand, it can be the actual probability with which the uncertainty parameter ξ takes on value in . On the other hand, P can simply describe the relative importance we assign to different uncertainty instances. We have the following definition: Definition 2.2 (Probability of Violation): Let x ∈ X be given. The probability of violation of x is defined as For example, if we assume a uniform probability density, then V(x) measures the ''volume of bad'' parameters ξ such that the constraint g(x, ξ ) ≤ 0 is violated. A solution x with small associated V(x) is feasible for most of the problem instances, i.e., it is approximately feasible for the worst-case problem. This concept of approximate feasibility has been introduced in the context of robust control in [15]. Any such solution is here named an '' -level'' solution: Definition 2.3 ( -Level Solution): Let ∈ (0, 1). We say that x ∈ X is an -level robustly feasible (or, more simply, Our ultimate goal is to devise an algorithm that returns a -level solution, where is any fixed small reliability level. The approach utilized in this paper uses a surrogate model called ''Scenario Design Problem''. By scenario it is here meant any possible realization or instance of the uncertainty parameter. In the ''scenario design'' we optimize the objective subject to a finite number of these randomly selected scenarios. Definition 2.4 (Scenario Design Problem): Assume that S independent identically distributed samples ξ 1 , . . . , ξ S are drawn according to probability P. A scenario design problem is given by the convex program The acronym SDP S refers to the fact that (2) is a convex program with S constraints. Here we assume the following technical condition on the scenario problem: The scenario problem SDP S is a standard convex optimization problem with a finite number of constraints S and, hence, its optimal solutionx S is (usually) efficiently computable by means of numerical algorithms [16].
The relationship between the number of sampled scenarios S and the probability of violation of the optimal solution to corresponding Scenario Design Problem V(x S ) was investigated in [17] and [18] -for a chosen we can always find S large enough such thatx S is -level feasible for the original problem (1) with arbitrarily high confidence. There is, however, no guarantee, that the resulting optimal objective value of (2) will be anywhere close to the true optimal value J * . The results derived in [18] show that the distribution function of V(x S ) is bounded by a beta distribution with parameters n x and S − n x + 1, and imposing that V(x S ) ≤ holds with high confidence implies that V(x S ) will be much less than in many cases, resulting in a conservative solution.
Next we introduce a concept that is crucial for the success of the Pooling part of the upcoming algorithm. Of the S generated scenarios, only some of these S will be ''bounding'' in the sense that they prevent the solution from ''falling'' to a lower objective value.
Definition 2.6 (Support Scenario): Scenario ξ i , i ∈ {1, . . . , S}, is a support scenario for the scenario problem SDP S if its removal changes the optimal solution of SDP S .
The following theorem, whose proof can be found in [18] or in a different form in [19], gives us the bound on the number of support scenarios:  The number of support scenarios for SDP S is at most n x , the size of x.
What is most important about this result is the fact that the number of support scenarios does not depend on the number of generated scenarios S. The first main contribution of this paper is an efficient way of solving (2), with the use of Theorem 2.7.

III. POOLING PART OF THE POOL & DISCARD ALGORITHM
The idea behind the Pooling part of the algorithm is the following: if one were to verbally describe the problem (2), the one word that came to our mind was ''long'', as there are much more constraints than decision variables. Moreover, the number of support constraints (or support scenarios), that the optimal solution of (2) depends upon is very small, when compared to the overall number of constraints (or scenarios).
The method consists of solving (2) by the following procedure. First, we start by completely neglecting the constraints in (2) that correspond to the different scenarios and solve this relaxed optimization problem. Then we find the most violated constraints (by computing the slacks), add them to the relaxed problem and find a new optimal solution.
The Pooling part can be summarized as follows: Step 0. Set I = ∅.
Step 1. Solve the following problem: and obtain a solutionx.
Step 2. Check feasibility of the solution by computing the slacks s i : Step 3. If max i∈{1,...,S} s i > 0, find the associated index of the maximum valueî = argmax i∈{1,...,S} s i , add it to the set I and return to Step 1. Otherwise, set x * =x, I * = I and terminate. It is important to remark that by the end of this procedure, we not only get the optimal solution of (2), but also an index set I that contains the support scenarios -this will be very significant for the success of the Discarding part of the P&D algorithm.
Another equally important remark concerns the efficient implementation of the algorithm. In Step 1, we are sequentially solving problems that are extremely similar, only differing in a single constraint. The use of warm-starts (when made possible by a proper choice of solution method) or even problem modification 1 (when supported by our choice of a solver) have immense effect on the efficiency of the Pooling part. For the implementation of the numerical examples that are investigated in this paper, we have chosen the JuMP package [20] for modeling optimization in the Julia language [21] and the CPLEX 12.7 solver [22]. This combination allowed us to use the algorithm to its full extent. 2 The machine, on which we conducted the numerical examples, was a PC with 3.6 GHz AMD Ryzen 5 2600X Six-Core CPU, 32 GB RAM, NVIDIA GeForce GTX 1050 Ti, running on 64-bit Windows 10.

A. NUMERICAL EXAMINATION -ASSET ALLOCATION PROBLEM
The first numerical example we chose to demonstrate the utility of the Pooling part of the P&D algorithm is the (by now, almost canonical) asset allocation problem [8]. Suppose we have n assets x 1 , . . . , x n that we want to invest in. The returns r 1 , . . . , r n of these assets are random variables. Our goal is to allocate our resources to these different assets, in order to maximize the quantile (often called the Value at Risk, or VaR) of the returns. This formulation neglects several of the important real-world issues -we do not allow short position, do not consider more than one trading period, etc. -the example is, above all else, intended to show the capabilities of the P&D algorithm. Our asset allocation problem can be summarized as follows: Our ability to solve (with no quotation marks) this problem depends heavily on the distribution of the returns r 1 , . . . , r n and the chosen quantile . Thanks to [23], we know that the feasible set of a scalar chance constraint is convex, provided that the vector (a T , b) T of the coefficients has symmetric logarithmically concave density and < 1/2. We will use this result and model the returns r as random variables that are independent and normally distributed (and, hence, have a symmetric logarithmically concave density). More precisely, the return r j has the following distribution i.e., the first return is ''deterministic'', with return r 1 = 1, and the nth return has mean µ n = 1.1 and standard deviation σ n = 0.1. Because of the chosen distribution of returns, the problem (5) can be transformed [8] into the following second order cone problem (SOCP, see [16]): where −1 (1 − ) is the 1 − quantile of the standard normal distribution. As an SOCP, this problem falls into the category of ''easy'' to solve (we can compute the optimal solution with little effort for large values of n -well into thousands) and as such provides the perfect ground for illustrating the capacities of the P&D algorithm. The scenario approach, works with a sample of S scenarios of the returns r i j , j = 1, . . . , n, i = 1, . . . , S. Using these scenarios, the sample counterpart to (5) has the following form: First of all, we will investigate on (7) the dependence of computation time of the Pooling part (CTPP) of the P&D algorithm for varying number of assets n and scenarios S. Additionally, we provide the computation time for the Pooling part without the use of warm-starts and problem modification (CnoWS) and the computation time for solving the problem (7) conventionally (CTC), i.e., passing it to the solver (CPLEX) with all the scenarios. The results of the computations are summarized in Table 1 and clearly demonstrate the effectiveness of the Pooling part of the P&D algorithm. As the number of scenarios grows, CTPP grows very slowly when compared to CTC, becoming over 20 times faster for the largest number of considered scenarios. The main factor in the effectiveness of the Pooling part is the low growth in the number of iterations needed to solve the problems with more scenarios -this should not be too surprising, since the number of support scenarios stays the same (for the same n). The variant without warm-start or problem modification CnoWS eventually (for high values of n) suffers from too big of an overhead when constructing the corresponding optimization problem, but can still outperform CTC in large number of instances.

IV. CONSTRAINT REMOVAL ALGORITHM
If all the S constraints are enforced, however, one cannot expect that good approximations of chance constrained solutions are obtained. To get a less conservative solution we use the framework introduced in [2] for relaxing problem (7). Their approach allows us to remove k constraints out of the S scenario constraints. A general removal procedure is formalized in the following definition: The sample-based optimization program where k constraints are removed as indicated by A is expressed as and its solution will be hereafter indicated as x * S,k . We introduce the following assumptions: This assumption requires that the algorithm A chooses constraints whose removal improves the solution by violating the removed constraints, and it rules out for example algorithms that remove inactive constraints only, or algorithms that remove constraints at random. Thus, this assumption is very natural and reflects the fact that we want to remove the constraints that improve the optimal objective value.
The next Theorem (proved in [2]) provides theoretical guarantees that V(x * S,k ) ≤ , i.e. that the optimal solution x * S,k of the optimization program SDP A S,k is feasible for the CCP .
then P S {V(x * S,k ) ≤ } ≥ 1 − β. The final result establishes that the objective value of CCP (whose optimal objective value will be denoted as J * ) can be approached at will, provided that sampled constraints are optimally removed. Let A opt be the optimal constraints removal algorithm which leads -among all possible eliminations of k constraints out of S -to the best possible improvement in the cost objective; further, let x * S,k,opt and J * S,k,opt be the corresponding optimal solution and objective value. We have the following theorem (again, proved in [2]).
Theorem 4.4 (Optimality): Let β ∈ (0, 1) be any small confidence parameter value, and let ν ∈ (0, ) be a performance degradation parameter value. If S and k are such that simultaneously hold with probability at least 1 − β.
One optimal way of removing constraints consists in discarding those constraints that lead to the largest possible improvement of the cost function. This approach is implemented by the following integer program, which has been described and investigated in [24], [25] and [26]: where M is a constant large enough so that, if z i = 1, then the constraint is satisfied for any candidate solution x. For k = 0, the formulations (2) and (11) are equivalent. By construction, problem (11) provides a framework for optimally selecting the constraints to be removed based on the inequality (10). However, solving (11) may be computationally challenging due to the increase in complexity from (2) to (11) that arises from the introduction of one binary variable per each of the S scenarios. In recent years, there have been developed strengthening procedures (see [27] and [28]) for some special structured problems, that significantly improve upon the formulation (11).

V. POOL & DISCARD ALGORITHM
The Discarding part of the algorithm consists of utilizing the index set I, finding the support scenarios among this set and finding the one scenario, whose removal decreases the optimal objective value the most -this is repeated k times, where k is either set a priori (by Theorem 4.3), or is terminated once an estimate of the probability of violation of obtained solution V(x) reaches certain threshold. This approach is almost identical to the one discussed in [29] (called greedy constraint removal), with the distinction that our algorithm utilizes the Pooling step and uses warm-starts (primarily utilizing I) throughout the iterations and as such can be rather effective (this will be demonstrated in the following sections). The P&D algorithm can be summarized as follows: Step 0. Solve the pooling part described above to obtain I * and x * . Set γ > 0, k > 0, I p = ∅. Repeat k times, or terminate once an estimate of V(x * ) reaches a threshold: Step 1. Find the set of support scenarios I r ⊂ I * -either by examining the slacks (s i > −γ ) or the associated dual variables (µ i > γ ).
Step 2. For each of the support scenarios i r ∈ I r , solve the following problem: using the Pooling part, warm-started by using I = I * \ {i r } and x = x * . Denote the solution to (12) VOLUME 8, 2020 as x * i r , its optimal objective function value v * i r and its final set of scenarios I * i r . Step 3. Find the index with the best optimal objective value: i * = argmin i r v * i r . Set x * = x * i * , I * = I * i * and add the corresponding scenario to the set of permanently discarded ones I p . The parameter γ can be, in theory, set to 0 -what discourages us from doing so are the implementation issues of numerical computing. When reporting the optimal dual variables µ the solvers rarely return exactly 0, more often, we get values ranging from 10 −8 to 10 −16 (the same goes for the slacks in the active constraints). If we did set γ to 0 we would (likely) have to consider all the scenarios as possible support scenarios and the execution of the algorithm would be significantly prolonged. Unless stated otherwise, the parameter γ was set to 10 −6 . It should be added, that Step 2. of the Discarding part can be fully parallelized to work more efficiently on multi-core machines or distributed computing environments.

A. NUMERICAL EXAMINATION -ASSET ALLOCATION PROBLEM CONTINUED
We will return to the same problem structure (7) again and examine the computational time for the whole P&D algorithm for varying number of variables and scenarios. In the Discarding part of the algorithm, we decided to discard k = S scenarios -note that this choice does not guarantee, that the resulting solution obtained by the P&D algorithm will be a -level feasible, not to mention having the objective value close to the optimal value objective J * .
To examine the effect of the warm-start in the discarding part (using the best solution x * and the index set I * from the previous iteration), we will first compare the computational times of the P&D algorithm, an algorithm that uses just the Pooling part without warm-starts (denoted as ''PnoD''), and an algorithm that uses neither Pooling nor Discarding (denoted as ''noPnoD'', which is essentially the one used in [29]), on a small-scale example (n = 20, = 0.01).
The comparison is summarized in Table 2 -the utilization of Pooling and the warm-starts in Discarding combined provide immense computational savings compared to the other two methods (while arriving at the exact same solution). To further compare the effectivity of the P&D algorithm, we set the parameters n, S, and k to the same values that can be found in [29] and compare the computation times directly (although they used different distributions for the asset returns, the problem structure is exactly the same). The results of the computation are reported in Table 3 -for n = 20, the results reported in [29] are comparable with the noPnoD variant of the algorithm, with slight improvement that is most likely caused by a more powerful machine and a newer version of the optimization solver. In the largest instance, the P&D algorithm was more that 200 times faster. For the n = 200, the authors in [29] used a random scenario removal strategy, instead of the greedy one (removing one of the support scenarios at random, instead of the one whose removal decreased the optimal objective value the most)this algorithm is O(n) times faster than the greedy one, but results in an inferior solution. In this setting, the P&D variant with randomized removal (denoted as P&D * ) was almost 500 times faster than the one in [29].
The real crux of the matter, however, is the following: ''How good a solution (in terms of -level feasibility and objective value) do we get by using the P&D algorithm?'' The remarkable thing about our optimal asset allocation problem is that for a chosen value of , we can get the optimal solution by solving the SOCP (6). Moreover, for every asset allocation x, we can find the corresponding quantile of the returns exactly. Or, alternatively, we can for a given value of the returns t and a given asset allocation x compute (again, exactly) the probability P{t ≤ n j=1 r j x j } (i.e., the smallest value of , for which our choice of x and t is feasible).
For the examination, we chose a problem with n = 30 assets and = 0.01. The optimal objective value (obtained by solving (6)) was 1.0309. The results are summarized in Table 4, Figure 1 and Figure 2. Using the formula for the needed number of scenarios from [18], with β = 10 −10 , we get that to obtain a feasible solution to this problem with high probability (1 − β), we need to solve (7) with at least S = 8,547 scenarios (without any discarding). The solution to this problem had the objective value 1.0179 (third column of Table 4), with the probability of violation 0.0029 (i.e. P{t ≤ n j=1 r j x j } = 0.0029) -i.e. we obtained a feasible solution, but with a rather poor objective value.
Afterwards, we ran the Discarding part of the algorithm, discarding S scenarios. The objective value improved to 1.0318 (fifth column of the table), but the corresponding probability of violation increased to 0.0138 -meaning that the combination of x and t (obtained after discarding) was no longer feasible. However, during the Discarding part of the algorithm we stored the particular solutions in each iteration. This allows us to find the last admissible (feasible) solution and find its corresponding objective and a number of scenarios that we discarded to get it -in this case the objective value was 1.0291 (seventh column of the table) with 31 discarded scenarios. An interesting thing to note is that even for 5,000 scenarios we still get a feasible solution (with probability of violation 0.0041) and can remove some scenarios, but for the lower numbers of scenarios even the ''robust'' solution is not feasible.
When we vary the number of scenarios several interesting phenomena appears. Firstly, when increasing the number of scenarios S we get a smaller value of the ''robust''   solution objective (the solution after the Pooling part) and smaller corresponding probability of violation (both of these are rather intuitive). Secondly, when we increase the number of scenarios, the probability of violation of the solution after discarding S scenarios approaches and the number of removed scenarios for an admissible solution gets closer to S . Thirdly, and most impressively, the admissible solution objective gets surprisingly close to the optimal value of (6).
Another feature of the P&D algorithm is that since we remove one scenario at a time, we can use the successive results to construct an approximation of the trade-off between reliability and optimal objective function value. This is best shown on Figure 1, where we can see the progression of the  P&D algorithm for different number of scenarios -each point corresponds to a solution with different number of removed scenarios (typically, more removed scenarios correspond to points more up and to the right). We included the optimal trade-off curve obtained by solving the SOCP (6) for different values of (called ''exact solution'' in the legend of Figure 2).
It must be emphasized that the P&D algorithm does not in any way incorporate any knowledge about the underlying distribution of the random variables. All it ''sees'' are the realizations in the form of individual scenarios.
The relationship between computational times, number of scenarios S, number of variables n and chosen probability of violation is further investigated in Tables 5, 6 and 7. From these results we can see that when we increase , the computation times increase linearly. The same cannot be said for when increasing S -the number of scenario removals and computational times of the pooling steps grow simultaneously, resulting in a superlinear increase computational time. Similar (and more impactful) behaviour can be observed in Table 7, where the number of variables n changes -this results in a larger number of possible support scenarios and larger solution times for the successive optimization problems.
When we increase the number of scenarios S, the resulting solutions (after discarding S scenarios) get very close to the true optimum J * , although it was not guaranteed by any theory. In similar fashion, the probability of violation of the solutions get very close to .

VI. NONLINEAR JOINT CHANCE CONSTRAINED EXAMPLE
In this section we investigate the performance of the algorithm on nonlinear example that appeared in the numerical sections of the state-of-the-art methods in [30] and [31]. Both of these methods are scenarios (or sample) based and use the indicator function approximation (although they approach it in different ways). In the method described in [30], the constraints need to be convex in x and the problem can be a joint chance constrained one. In the method described in [31], the constraints do not have to be convex, but must be continuously differentiable in x and the authors deal with a single chance constraint only. The problem both papers have chosen for the numerical examination is the following one: where ξ ij , i = 1, . . . , m and j = 1, . . . , n are independent and identically distributed standard normal random variables, b ∈ . In the case of [31], m = 1 (a single chance constraint).
Optimal solution x * of the problem (13), derived in [30], is: where F −1 χ 2 n is the inverse chi-squared distribution function with n degrees of freedom.
Because of the nature of the problem (quadratic and convex), we were able to use the CPLEX solver, and utilize the problem modification feature again. We start the numerical examination with the same setting as [31]: n = 10, m = 1, b = 10, = 0.05. The parameter γ that controls the selection of scenarios to discard, was set to 10 −3 . Using the formula (14), the optimal objective value of this problem is −7.390. We generate a number of scenarios S (the values were log-spaced between 10 2 and 10 4 ) and set the P&D algorithm to discard S of them. After that we estimate the reliability (1 − ) of the obtained solution using 10 5 new scenarios. The results of the computations are summarized in Table 8. Unsurprisingly, the more scenarios are taken into account, the better (closer to the theoretical optimum) the result. The computational time is quite good, considering [31] report around 500 s as the computational time for their algorithm (that uses 500 scenarios and reports the optimal objective value -7.627) and the big-M mixed-integer formulation does not converge in an hour [31] (again, using ''just'' 500 scenarios).
The second setting we investigate is from [30]: n = 10, m = 10, b = 100, = 0.1. The optimal objective value, using (14), is −20.82. Note that in this setting we are dealing with a ''proper'' joint chance constraint problem, with nonlinear (but convex) constraint functions. In the Pooling part of the algorithm, when we find a scenario with a violated constraint, we have a choice of either adding all of the m constraints that correspond to this scenario to the problem, or to add only the violated ones. In our implementation we chose the former, since it corresponds more closely to the description of the P&D algorithm we gave in the previous sections, although additional computational savings could be gained by properly implementing the latter approach. The number of scenarios used in the examination ranged between 10 2 and 10 4 and the number of scenarios to discard was set to S . The results of the computations are listed in Table 9. To compare the results with the ones achieved in [30], where they used 10,000 scenarios for the computations -the numbers the authors report are a bit vague: ''Our algorithm typically requires less than 10 iterations to converge to the optimal value, and each iteration approximately takes 6 s on average.'' The main objections being that their algorithm was presented on two problems with different dimensions (the one presented here and a smaller one), and that, at least judging from the figures (as there is no other way to find the value), their ''optimal solution'' was around −20.4, which is rather far from the real one. It is important to emphasize again that the P&D algorithm does not produce just one solution -as a sort of a by-product it generates a sequence of decisions, that are ''optimal'' with respect to an increasing number of discarded scenarios. When we estimate the reliability of these solutions, we get an approximation of the trade-off between the reliability level (1-) and optimal objective value. Naturally, this approximation gets better as we increase the number of scenarios. The approximation of the trade-off for the setting described above is depicted in Figure 3 -it shows just one run of the P&D algorithm for different number of scenarios and the optimal values computed using the formula (14). The reliability of the solution is estimated using 10 5 different scenarios. If we stopped the algorithm with 200 scenarios once the estimate of the reliability of the solution gets lower than the desired level 1 − and use the previous value, we would discard only 12 scenarios with the objective value −20.47 and estimated reliability 0.9143 -this takes around 6 s. Note that the robust solution for this problem, i.e. the solution for = 0, is clearly 0, since each ξ 2 ij can attain any nonnegative value.

VII. CONCLUSION
The main advantage of the P&D algorithm lies in the exploitation of the structure of the scenario design problem, which is done on two levels: • The Pooling part of the P&D algorithm utilizes the fact that the number of support scenarios is usually very small compared to the number of sampled scenarios. By iteratively solving much smaller problems we can get the solution faster and need less memory, than we would need to solve the scenario design problem with all scenarios at once (compare the columns CTC and CTPP in Table 1).
• The Discarding part of the algorithm fully utilizes the set I that contains the current support scenarios, to find the ones that will be discarded (either be the greedy or the randomized algorithm). Since it only solves comparatively much smaller problems (and, each scenario removal should terminate in just a few iterations), it brings additional computational savings (compare the columns PnoD and P&D in Table 2). The combined effect of the two parts of the algorithm is best seen in the difference between columns noPnoD and P&D in Table 2. The numerical examinations show that P&D algorithm provides a powerful framework for handling certain types of chance constrained optimization problems. When compared with conventional approaches [29] on a linear example (see Table 3), it was several hundred times more efficient in the largest instances. On the nonlinear examples, it was on par with the state-of-the-art methods [31] and [30].
Further investigation are possible -in the Pooling part of the algorithm for joint chance constrained problems, the choice between including all constraints for a violated scenario, or just the ones that are violated, could bring additional computational savings. Also, the use of P&D (or just the use of Pooling) could be applied in other classes of convex optimization problems (e.i., semi-definite problems in control) that need to be set in chance constrained (or robust) setting and require the consideration of large number of scenarios.
JAKUB KŮDELA received the M.S. degree in mathematical engineering from the Brno University of Technology, in 2014, and the Ph.D. degree in applied mathematics from the Brno University of Technology in 2019.
Since 2018, he has been a Research Assistant with the Institute of Automation and Computer Science, Brno University of Technology. His research interests include the development of computational methods for various optimization problems and engineering applications.
PAVEL POPELA received the Ph.D. degree in econometrics from Charles University Prague, in 1998. Since 1986, he has been as a Senior Lecturer and a Researcher with the Brno University of Technology. His research interests include stochastic programming models and methods. VOLUME 8, 2020