Runtime Analysis of Single- and Multi-Objective Evolutionary Algorithms for Chance Constrained Optimization Problems with Normally Distributed Random Variables

1 Introduction

Many real-world optimization problems involve solving optimization problems that contain stochastic components [BenTal et al., 2009]. Chance constraints [Charnes and Cooper, 1959] allow to limit the probability of violating a constraint involving stochastic components. In contrast to limiting themselves to ruling out constraint violations completely, this allows to deal with crucial constraints in a way that allows to ensure meeting the constraints with high confidence (usually determined by a confidence level α) while still maintaining solutions of high quality with respect to the given objective function.

Evolutionary algorithms have successfully been applied chance constrained optimization problems [Poojari and Varghese, 2008; Liu et al., 2013]. Recent studies investigated the classical knapsack problem in static [Xie et al., 2019; Xie et al., 2020; Neumann et al., 2022] and dynamic settings [Assimi et al., 2020] as well as complex stockpile blending problems [Xie et al., 2021a] and the optimization of submodular functions [Neumann and Neumann, 2020]. Theoretical analyses for submodular problems with chance constraints, where each stochastic component is uniformly distributed and has the same amount of uncertainty, have shown that greedy algorithms and evolutionary Pareto optimization approaches only lose a small amount in terms of approximation quality when comparing against the corresponding deterministic problems [Doerr et al., 2020; Neumann and Neumann, 2020] and that evolutionary algorithms significantly outperform the greedy approaches in practice. Other recent theoretical runtime analyses of evolutionary algorithms have produced initial results for restricted classes of instances of the knapsack problem where the weights are chosen randomly [Neumann and Sutton, 2019; Xie et al., 2021b].

For our theoretical investigations, we use runtime analysis which is a major theoretical tool for analyzing evolutionary algorithms in discrete search spaces [Neumann and Witt, 2010; Jansen, 2013; Doerr and Neumann, 2020]. In order to understand the working behaviour of evolutionary algorithms on broader classes of problems with chance constraints, we consider the optimization of linear functions with respect to chance constraints where the stochastic components are independent and each weight wi is chosen according to a Normal distribution N(µi, σ2 i ). This allows to reformulate the problem by a deterministic equivalent non-linear formulation involving a linear combination of the expected expected value and the standard deviation of a given solution.

We investigate how evolutionary algorithms can deal with chance constrained problems where the stochastic elements follow a Normal distribution. We first analyze the classical (1+1) EA on this problem formulation and we show that imposing a simple cardinality constraint for a simplified class of instances leads to local optima and exponential lower bounds for the (1+1) EA. In order to deal with the issue of the (1+1) EA not being able to handle even simple constraints due to the non-linearity of the objective functions, we introduce a Pareto optimization approach for the chance constrained optimization problems under investigation. So far, Pareto optimization approaches that achieved provably good solution provided a trade-off with respect to the original objective functions and given constraints. In contrast to this, our approach trades off the different components determining the uncertainty of solutions, namely the expected value and variance of a solution. A crucial property of our reformulation is that the extreme points of the Pareto front provide optimal solutions for any linear combination of the expected value and the standard deviation and solves the original chance constrained problem for any confidence level α ≥ 1/2. These insights mean that the users of the evolutionary multi-objective algorithm does not need to know the desired confidence level in advance, but can pick from a set of trade-offs with respect to the expected value and variance for all possible confidence levels. We show that this approach can also be applied to the chance constrained minimum spanning tree problem where each edge cost is chosen independently according to its own Normal distribution. In terms of algorithms, we analyze the well-known GSEMO [Giel, 2003] which has been frequently applied in the context of Pareto optimization [Neumann and Wegener, 2006; Kratsch and Neumann, 2013; Friedrich and Neumann, 2015; Zhou et al., 2019] and show that it computes such an optimal set of solutions for any confidence level of α ≥ 1/2 in expected polynomial time if the population size stays polynomial with respect to the given inputs. Finally, we experimentally compare the (1+1) EA and GSEMO on different stochastic instances of the NP-hard minimum weight dominating set problem and show that GSEMO significantly outperforms the (1+1) EA for almost all problem instances.

2 Chance Constrained Optimization Problems

Our basic chance-constrained setting is given as follows. Given a set of n items E = {e1, . . . , en} with weights wi, 1 ≤ i ≤ n, we want to solve

where w(x) = ∑ni=1 wixi, x ∈ {0, 1}n , and α ∈ [0, 1]. Throughout this paper, we assume that the weights are independent and each wi is distributed according to a Normal distribution N(µi, σ2 i ), 1 ≤ i ≤ n, where µi ≥ 1 and σi ≥ 1, 1 ≤ i ≤ n. We denote by µmax = max1≤i≤n µi and vmax = max1≤i≤n σ2 i the maximal expected value and maximal variance, respectively. According to [Ishii et al., 1981], the problem given in Equation 1 is in this case equivalent to minimizing

where Kα is the α-fractile point of the standard Normal distribution. Throughout this paper, we assume α ∈[1/2, 1[ as we are interested in solutions of high confidence. Note that there is no finite value of Kα for α = 1 due to the infinite tail of the Normal distribution. Our range of α implies Kα ≥ 0.

We carry out our investigations where there are additional constraints. First, we consider the additional constraint

Algorithm 1: (1+1) EA for minimization

1 Choose x ∈ {0, 1}n uniformly at random; 2 repeat 3 Create y by flipping each bit xi of x with probability 1n ; 4 if f(y) ≤ f(x) then 5 x ← y;

• 6 until stop;

|x|1 ≥ k, which requires that at least k items are chosen in each feasible solution. Furthermore, we consider the formulation of the stochastic minimum spanning tree problem given in [Ishii et al., 1981]. Given an undirected connected weighted graph G = (V, E), n = |V | and m = |E| with randomweights w(ei), ei ∈ E. The search space is {0, 1}m. For a search point x ∈ {0, 1}m, we have w(x) = ∑mi=1 w(ei)xi as the weight of a solution x. We investigate the problem given in Equation 1 and require for a solution x to be feasible that x encodes a connected graph. We do not require a solution to be a spanning tree in order to be feasible as removing an edge from a cycle in a connected graph automatically improves the solution quality and is being taken care of by the multi-objective algorithms we analyze in this paper. Note that the only difference compared to the previous setting involving the uniform constraint is the requirement that a feasible solution has to be a connected graph.

3 Analysis of (1+1) EA

The (1+1) EA (Algorithm 1) is a simple evolutionary algorithm using independent bit flips and elitist selection. It is very well studied in the theory of evolutionary computation [Doerr, 2020] and serves as a stepping stone towards the analysis of more complicated evolutionary algorithms. As common, in the area of runtime analysis, we measure the runtime of the (1+1) EA by the number of iterations of the repeat loop. The optimization time refers to the number of fitness evaluations until an optimal solution has been obtained for the first time, and the expected optimization time refers to the expectation of this value.

3.1 Lower Bound for (1+1) EA and Uniform Constraint

We consider the (1+1) EA for the problem stated in Equation 1 with an additional uniform constraint that requires that each feasible solution contains at least k elements, i.e. |x|1 ≥ k holds. We show that the (1+1) EA has an exponential optimization time on an even very restrictive class of instances involving only two different weight distributions.

We use the following fitness function, which should be minimized in the (1+1) EA:

where L = (1 + ∑ni=1 µi + Kα(∑ni=1 σ2 i )1/2 ). This gives a large penalty to each unit of constraint violation. It implies that any feasible solution is better than any infeasible solution and that f(x) > f(y) holds if both x and y are infeasible and |x|1 > |y|1. Furthermore, the fitness value of an infeasible solution only depends on the number of its elements.

We now show a lower bound on the optimization time of the (1+1) EA for a specific instance class I containing only two types of elements. Type a elements have weights chosen according to N(n2 + δ, 1) and type b elements have weights chosen according to N(n2 , 2). We set δ = 1 2√k·1.48 . The instance I has exactly n/2 elements of type a and n/2 elements of type b. We consider Kα = 1 which matches α ≈ 0.84134, and set k = 0.51n. Using the fitness function f, we have the additional property for two feasible solutions x and y that f(x) < f(y) if |x|1 < |y|1 due to an expected weight of at least n2 for any additional element in a feasible solution. This also implies that an optimal solution has to consist of exactly k elements. The quality of a solution with k elements only depends on the number of type a and type b elements it contains. An optimal solution includes n/2 elements of type a whereas a locally optimal solution includes n/2 elements of type b. Note that an optimal solution has minimal variance among all feasible solutions. The (1+1) EA produces with high probability the locally optimal solution before the global optimum which leads to the following result.

Theorem 1. The optimization time of the (1+1) EA on the instance I using the fitness function f is eΩ(n) with probability 1 - e-Ω(n1/4 ) .

A similar lower bound for a specific class of instances of the chance constrained minimum spanning tree problem having two types of edge weights can be obtained.

We now show the asymptotic behaviour of the (1+1) EA on instance I through an experimental study. Table 1 shows for n ∈ {100, 200, 500, 1000, 1500, 2000} the number of times out of 30 runs the globally optimal solution has been obtained before the locally optimal one. Note that it takes exponential time to escape a locally optimal solution. It can be observed that the fraction of successful runs obtaining the global optimum clearly decreases with n. For n = 2000 no globally optimal solution is obtained within 30 runs.

4 Multi-Objective Evolutionary Algorithm

We now introduce bi-objective formulations of the chance constrained problems with uniform and spanning tree constraints. We use a Pareto Optimisation approach for this which computes trade-offs with respect to the expected weight µ and variance v. We say that a solution z dominates a solution x (denoted as z ≼ x) iff µ(z) ≤ µ(x) and

Algorithm 2: Global SEMO

1 Choose x ∈ {0, 1}n uniformly at random; 2 P ← {x}; 3 repeat 4 Choose x ∈ P uniformly at random; 5 Create y by flipping each bit xi of x with probability 1n ; 6 if ∄w ∈ P : w ≺ y then 7 S ← (P ∪ {y})\{z ∈ P | y ⪯ z};

8 until stop;

v(z) ≤ v(x). We say that z strongly dominates x (donated as z ≺ x) iff z ⪯ x and µ(z) < µ(x) or v(z) < v(x).

We investigate the algorithm GSEMO [Laumanns et al., 2004; Giel, 2003] shown in Algorithm 2, which has been frequently used in theoretical studies of Pareto optimization. It starts with a solution chosen uniformly at random and keeps at each time step a set of non dominated solutions found so far. In addition to being able to achieve strong theoretical guarantees [Friedrich et al., 2010; Friedrich and Neumann, 2015; Zhou et al., 2019], GSEMO using different types of multi-objective formulation has shown strong performance in practice [Qian et al., 2015; Qian et al., 2017; Roostapour et al., 2022]. We study the multi-objective evolutionary algorithms in terms of the expected time (measured in terms of iterations of the algorithm) until they have produced a population which contains an optimal solution for each α ∈ [1/2, 1[.

4.1 UniformConstraints

For the case of the uniform constraint |x|1 ≥ k, we consider the objective function f(x) = (µ(x), v(x)) where

Note that it gives the expected value and variance for any feasible solution, and a large penalty for any unit of constraint violation in each objective function if a solution is infeasible. This implies that the objective value of an infeasible solution is always worse than the value of a feasible solution.

As we have µi ≥ 1 and σi ≥ 1, 1 ≤ i ≤ n, each Pareto optimal solution contains exactly k elements. This is due to the fact that we can remove from a solution x with |x|1 > k any element to obtain a solution y with µ(y) < µ(x) and v(y) < v(x). We will minimize fλ(x) = λµ(x)+(1-λ)v(x) by selecting minimal elements with respect to fλ(ei) = λµi + (1 - λ)σ2 i , 0 < λ < 1. For the special cases λ = 0 and λ = 1, we minimize fλ by minimizing f0(x) = (v(x), µ(x)) and f1(x) = (µ(x), v(x)) with respect to the lexicographic order. Note, that we are using fλ both for the evaluation of a search point x as well as the evaluation of an element ei. For each fixed λ ∈ [0, 1], an optimal solution for fλ can be obtained by sorting the items increasing order of fλ and selecting the first k of them. For a given set X of such points we denote by X∗ λ ⊆ X the set of minimal elements with respect to fλ. Note that all points in the sets X∗ λ , 0 ≤ λ ≤ 1, are not strongly dominated in X and therefore constitute Pareto optimal points when only considering the set of search points in X.

Definition 1 (Extreme point of set X). For a given set X, we call f(x) = (µ(x), v(x)) an extreme point of X if there is a λ ∈ [0, 1] such that x ∈ X∗ λ and v(x) = maxy∈X∗ λ v(x).

We denote by f(X) the set of objective vectors corresponding to a set X ⊆ 2E , and by f(2E ) the set of all objective vectors of the considered search space 2E . The extreme points of 2E are given by the extreme points of the convex hull of f(2E ). A crucial property of the extreme points is that they contain all objective vectors that are optimal for any λ ∈ [0, 1]. Hence, if there is an optimal solution that can be obtained by minimizing fλ for a (potentially unknown) value of λ, then such a solution is contained in the set of search points corresponding to the extreme points of 2E .

In the following, we relate an optimal solution of

subject to |x|1 ≥ k.

to an optimal solution of

for a given parameter R ≥ 0 that determines the weightening of µ(x) and v(x). Note that gR is a linear combination of the expected value and the variance and optimizing gR is equivalent to optimizing fλ for λ = R/(R+Kα) as we have gR(x) = (R+Kα) · fλ(x) in this case. We use gR to show that there is a weightening that leads to an optimal solution for g following the notation given in [Ishii et al., 1981], but will work with the normalized weightening of λ when analyzing our multi-objective approach.

Let x∗ be an (unknown) optimal solution for g and let D(x∗ ) = (∑ni=1 σ2 i x∗i )1/2 be its standard deviation. Lemma 1 follows directly from the proof of Theorems 1-3 in [Ishii et al., 1981] where it has been shown to hold for the constraint where a feasible solution has to be a spanning tree. However, the proof only uses the standard deviation of an optimal solution and relates this to the weightening of the expected value and the variance. It therefore holds for the whole search space independently of the constraint that is imposed on it. Therefore, it also holds for the uniform constraint where we require |x|1 ≥ k.

Lemma 1 (follows from Theorems 1-3 in [Ishii et al., 1981]). An optimal solution for g2D(X∗ ) is also optimal for g.

Based on Lemma 1, an optimal solution for fλ, where λ = 2D(X∗ )/(2D(X∗ ) + Kα), is also optimal for g. As we are dealing with a uniform constraint, an optimal solution for fλ can be obtained by greedily selecting elements according to fλ until k elements have been included. The extreme points of the convex hull allow to cover all values of λ where optimal solutions differ as they constitute the values of λ where the optimal greedy solution might change and we bound the number of such extreme points in the following.

In order to identify the extreme points of the Pareto front, we observe that the order of two elements ei and ej with respect to a greedy approach selecting always a minimal element with respect to fλ can only change for one fixed value

of λ. We define λi,j = σ2 j -σ2 i (µi-µj )+(σ2 j -σ2 i ) for the pair of items ei and ej where σ2 i < σ2 j and µi > µj holds, 1 ≤ i < j ≤ n.

Consider the set Λ = {λ0, λ1, . . . , λℓ, λℓ+1} where λ1, . . . , λℓ are the values λi,j in increasing order and λ0 = 0 and λℓ+1 = 1. The key observation is that computing Pareto optimal solutions that are optimal solutions for fλ and every λ ∈ Λ gives the extreme points of the problem.

Lemma 2. Each extreme point of the multi-objective formulation is Pareto optimal and optimal with respect to fλ for at least one λ ∈ Λ. The number of extreme points is at most n(n - 1)/2 + 2 ≤ n2 .

In the following, we assume that vmax ≤ µmax holds. Otherwise, the bound can be tightened by replacing vmax by µmax. The following lemma gives an upper on the expected time until GSEMO has obtained Pareto optimal solution of minimal variance. Note that this solution is optimal for f0. We denote by Pmax the maximum population size that GSEMO encounters during the run of the algorithm.

Lemma 3. The expected time until GSEMO has included a Pareto optimal search point of minimal variance in the population is O(Pmaxn2 (log n + log vmax)).

Lemma 3 directly gives an upper bound for GSEMO on the worst case instance I for which the (1+1) EA has an exponential optimization time. The feasible solution of minimal variance is optimal for instance I. Furthermore, we have vmax ≤ 2n and the variance can only take on O(n) different values. which implies that Pmax = O(n) holds for instance I. Therefore, we get the following result for GSEMO on the worst case instance I of the (1+1) EA.

Theorem 2. The expected time until GSEMO has obtained an optimal solution for the instance I is O(n3 logn).

We now present the result for the general class with a uniform constraint. Based on a Pareto optimal solution having the minimal variance, GSEMO can construct all other extreme points and we obtain the following result.

Theorem 3. Considering the chance constrained problem with a uniform constraint, the expected time until GSEMO has computed a population which includes an optimal solution for any choice of α ∈ [1/2, 1[ is O(Pmaxn2 ℓ(logn + log vmax)).

Proof. We assume that we have already included a Pareto optimal solution of minimal variance vmin into the population. Let

and

be the maximal and minimal variance of any optimal solution for the linear weightening fλ.

Note that we have vmax λ = vmin λ for λ = 0 as the Pareto optimal objective vector of minimal variance is unique. Hence, the Pareto optimal solution of minimal variance vmin = vmax 0 is a solution of maximal variance for λ = 0. Consider λi, 0 ≤ i ≤ ℓ. We have vmax λi = vmin λi+1 as the extremal point that is optimal for fλi and fλi+1 has the largest variance for fλi and the smallest variance for fλi+1 among the corresponding sets of optimal solutions.

Assume that we have already included into the population a search point x that is minimal with respect to fλi and has maximal variance vλmax i among all these solutions. The solution x is also optimal with respect to fλi+1 and we have vmax λi = vmin λi+1 . We let r be the number of elements contained in x but not contained in the optimal solution y for fλi+1 that has variance vλmax i+1 . As both solutions contain k elements, differ by r elements and y has maximal variance with respect to fλi+1 , there are r2 2-bit flips that bring down the distance d(x) = vmax λi+1 - v(x) ≤ vmax λi+1 - vmax λi . Using the multiplicative drift theorem [Doerr et al., 2012] where we always choose the solution that is optimal with respect fλi+1 and has the largest variance, the expected time to reach such as solution of variance vλmax i+1 is O(Pmaxn2 log(vmax λi+1 - vmax λi )) = O(Pmaxn2 (logn + log vmax)). Summing up over the different values of i, we get O(Pmaxn2 ℓ(log n + log vmax)) as an upper bound on the expected time to generate all extreme points.

4.2 Chance ConstrainedMinimum Spanning Trees

We now extend the previous results to the chance constrained minimum spanning tree problem where edge weights are independent and chosen according to a Normal distribution. Note that using the expected weight and the variance of a solution as objectives results in a bi-objective minimum spanning tree problems for which a runtime analysis of GSEMO has been provided in [Neumann, 2007].

Let c(x) be the number of connected components of the solution x. We consider the bi-objective formulation for the multi-objective minimum spanning tree problem given in [Neumann, 2007]. Let wub = n2 · max{µmax, vmax}. The fitness of a search point x is given as f(x) = (µ(x), v(x)) where µ(x) = (c(x) - 1) · wub + ∑mi=1 µixi and v(x) = (c(x) - 1) · wub + ∑mi=1 σ2 i xi. It gives a large penalty for each additional connected component.

We transfer the results for the multi-objective setting under the uniform constraint to the setting where a feasible solution has to be a spanning tree. The crucial observation from [Neumann, 2007] to obtain the extreme points is that edge exchanges resulting in new spanning trees allow to construct solutions on the linear segments between two consecutive extreme points in the same way as in the case of the uniform constraint. Let ℓ ≤ m(m - 1)/2 be the pairs of edges ei and ej with σ2 i < σ2 j and µi > µj. Similar to Lemma 2 and using the arguments in [Ishii et al., 1981], the number of extreme points is at most ℓ + 2 ≤ m2 as an optimal solution for fλ = λµ(x) + (1 - λ)v(x) can be obtained by Kruskal's greedy algorithm. We replace the expected time of O(n2 ) for an items exchange in the case of the uniform constraint with the expected waiting time of O(m2 ) for a specific edge exchange in the case of the multi-objective spanning tree formulation and get the following results.

Theorem 4. Considering the chance constrained minimum spanning tree problem, the expected time until GSEMO has computed a population which includes an optimal solution for any choice of α ∈ [1/2, 1[ is O(Pmaxm2 ℓ(logn + log vmax)).

5 Experimental Investigations

We now present experimental results for the chance constrained version of a classical NP-hard optimization problem. We consider the minimum weight dominating set problem. Given a graph G = (V,E) with weights on the nodes, the goal is to compute a set of nodes D of minimal weight such that each node of the graph is dominated by D, i.e. either contained in D or adjacent to a node in D. Let n = |V | be the number of nodes in the given graph G = (V, E). To generate the benchmarks, we assign each node u ∈ V a Normal distribution N(µ(u), v(u)) with expected weight µ(u) and a variance v(u). We consider for each graph values of α = 1 - β where β ∈ {10-2 , 10-4 , 10-6 , 10-8 , 10-10 }. The objective function is given by Equation 2 where we consider the chosen nodes together with their weight probability distributions. A solution is feasible if it is a dominating set. Therefore, we use the fitness functions from Section 3.1 for the (1+1) EA and Section 4.1 for GSEMO and have a large penalty term as before for each node in V that is not dominated in the given search point x.

We investigate the graphs cfat200-1, cfat200-2, ca-GrQC and Erdoes992, which are sparse graphs chosen from the network repository [Rossi and Ahmed, 2015]. Each expected weight µ(u) is an integer chosen independently and uniformly at random in {n2 , . . . , 2n2 } and the variance v(u) is an integer chosen independently and uniformly at random in {n3 , . . . , 2n3 } for cfat200-1 and cfat200-2 and an integer chosen uniformly at random in {n, . . . , n2 } for ca-GrQC and Erdoes992. For the graphs cfat200-1 and cfat200-2 (out of the DIMACS benchmark set), which consist of 200 nodes each, we give each algorithm a budget of 3n3 log n fitness evaluations. This choice for a relatively small number of nodes is based on common results in the area of runtime analysis. It assumes that the population size is of the order n and that 2-bit flips (costing a factor of n2 ) involving a coupon collector effect (factor logn) constitute the essential part of the optimization process. The constant "3" is a bit larger than the usual factor "e" (Eulerian number) often seen in runtime analyses and gives the algorithms a slightly larger budget which would reduce failure probabilities during the run of an algorithm. For the graphs ca-GrQc and Erdos992 (out of the collaboration network benchmark set), which consist of 4158 and 6100 nodes, respectively, we allocate a budget of 10 million fitness evaluations to each run. The goal is here to investigate the quality of the results obtained within a reasonable time budget. GSEMO spends its whole fitness budget on a single run whereas for the (1+1) EA the budget is equally divided among the 5 runs for the different 5 values of α. For each graph and α combination we obtain 30 results as described before. Note that the although GSEMO obtains in one run results for each value of α, the results for a fixed graph and α combination are independent of reach other. We report on the mean and standard deviation of the result obtained for (1+1) EA and GSEMO and use the Mann-Whitney test to compute the p-value in order to examine statistical significance. We call a result statistically significant if the p-value is at most 0.05.

The results are shown in Table 2. It can be observed that GSEMO outperforms the (1+1) EA for 18 out of the 20 settings in terms of the mean value that is achieved in 30 runs. All results are statistically significant. Overall, this shows a clear advantage of using the multi-objective model presented in this paper. We finally report on the maximum population size Pmax during the runs as this is an important parameter in the runtime analysis we carried out. Our experiments show that the maximum population is in general small and does not grow exponentially with the number of nodes of the given graph. The maximum population size of GSEMO during any of the 30 runs has been 127 for cfat200-1, 34 for cfat200-2, 15 for ca-GrQc, and 7 for Erdos992. This shows the the maximum population size stays moderate for the problem instances we investigated experimentally.

6 Conclusions

With this paper, we provided the first analysis of evolutionary algorithms for chance constrained combinatorial optimization problems with Normally distributed variables. For the case of uniform constraints we have shown that there are simple instances where the (1+1) EA has an exponential optimization time. Based on these insights we a multi-objective formulation which allows Pareto optimization approaches to compute a set of solutions containing for every possible confidence of α an optimal solution. Finally, we showed the effectiveness of our multi-objective approach by investigating the chance constrained setting of the minimum weight dominating set problem.

Acknowledgments

This work has been supported by the Australian Research Council (ARC) through grant FT200100536 and by the Independent Research Fund Denmark through grant DFF-FNU 8021-00260B.