Improved Integer Programming Approaches for Chance-Constrained Stochastic Programming

1 Introduction

Decision making is one of the important areas in artificial intelligence. Not to mention the frame problem [McCarthy and Hayes, 1969], when we make decisions, it is unrealistic to assume that all of the necessary information is available (e.g., [Bonet andGeffner, 2011]). Therefore, it is important to model the decision making problem as an appropriate model, which depends on the available information. For example, we can model a problem with some random variables as a stochastic programming problem [Heyman and Sobel, 2003] if we know the distributions of the random variables. For another example, we can model a problem with some uncertain values as a robust optimization problem [Ben-Tal et al., 2009] if we know the maximum and minimum values of the uncertain values.

A Chance-Constrained Stochastic Programming (CCSP) problem is an important stochastic programming problem having many applications, where an objective function is minimized under probabilistic constraints. It can be formulated as

where c is a vector, ˜ A and ˜b are a random matrix and a vector, respectively, ϵ is a constant parameter between 0 and 1, and X is the solution space of the decision variable x. The objective of the CCSP is to find x ∈ X that minimizes the inner product cx under the probabilistic constraint that the relation ˜ Ax ≥ ˜b holds with probability at least 1 - ϵ. Applications of CCSP problems include production planning [Lejeune and Ruszczy´nski, 2007] and the probabilistic lot-sizing problem [Beraldi and Ruszczy´nski, 2002].

The studies for CCSP have a long history. Charnes, Cooper, and Symonds [1958] first studied the CCSP with disjoint probabilistic constraints, where the probability of violating each constraint is separately evaluated. The CCSP with joint probabilistic constraints was first studied by Miller and Wagner [1965], but the random variables in their formulation are independent. The CCSP with joint probabilistic constraints with dependent random variables was first introduced in Pr´ekopa [1973].

In this paper, we consider a special case of the CCSP, where (i) the random variables appear only in the righthand side of the constraints and (ii) the random variable ˜b has a discrete distribution having a finite support. Condition (i) means that the probabilistic constraint can be written as Pr(Ax ≥ ˜b) ≥ 1 - ϵ, where A is a constant matrix. Condition (ii) means that each random variable in ˜b has a set of discrete realizations with probabilities. If we allow approximation, Condition (ii) can be avoided by converting each random variable with an infinite or continuous support into a discrete random variable with a finite number of realizations using sampling methods. Since there are many applications of the special case of the CCSP (such as the unit commitment problem), there exist many papers on this constrained CCSP such as [Luedtke et al., 2010; K¨uc¸ ¨ukyavuz, 2012].

A standard approach to obtain the optimal solution for the CCSP (introduced in [Ruszczy´nski, 2002]) is to formulate it as a Mixed Integer Programming (MIP), but this requires scenario enumeration. This means that we need to enumerate all of the possible combinations of the realizations of the random variables (i.e., scenarios) to formulate the MIP. If we enumerate all of the scenarios, then the number of scenarios is exponentially large with regard to the number of random variables. Since the size of the MIP problem is proportional to the number of scenarios, the resulting MIP formulation becomes huge. In practice, we cannot exactly solve such a large MIP formulation.

In this paper, we show how to reduce the number of required scenarios by using a vector domination relation without loss of accuracy. While the standard approach includes all of the scenarios S in the MIP formulation, our new idea is to construct an MIP formulation that includes only a small subset Sϵ ⊆ S. Although the number of scenarios |Sϵ| in the worst case is still exponential in the number of random variables, the number of scenarios is significantly reduced in practice, which allows us to exactly solve a larger instance of the CCSP.

In addition, we propose an approach to obtain an approximate solution for huge CCSP problems, by using Boole's inequality (also known as the union bound) to construct a compact formulation of the CCSP, avoiding any scenario generation. Similar to the scenario reduction, our new approach is to identify an approximate subset S ′ϵ ⊆ Sϵ such that S ′ϵ and Sϵ are similar, using Boole's inequality. Then we formulate a new MIP problem using S ′ϵ . This approach can be seen as decomposing the joint probabilistic constraints in the original CCSP into a problem with individual chance constraints.

We also show that the number of binary variables in the compact CCSP formulation can be further reduced without loss of accuracy. Using technique similar to the ones in [Ibaraki, 1976; Vielma et al., 2010], we reformulate the compact CCSP problem by reducing the number of binary variables logarithmically.

2 The Unit Commitment Problem

The unit commitment problem is one of the important applications of a special case of the CCSP in which only the right-hand side vector is random with a discrete distribution having a finite support. In electric power utility companies, they want to turn off unneeded generators to save cost. However, the generators cannot be switched quickly between the on and off states, and the switching is expensive. Therefore typical electric utilities use predictions of the electricity demands and create a schedule for turning the generators on and off. Ignoring some of the complex constraints arising in real schedules of generators, we can formulate naturally the most essential part of the unit commitment problem as this CCSP

where the constant gi is the maximum power generator i can generate, the constant ci is the cost for using generator i, and the random variable ˜ dt is the power demand at time t. The decision variable xit is equal to 1 if the generator i is active at time t (that is, it is generating electricity), and otherwise inactive (that is, the generator is turned off). The objective is to minimize the total cost for using these power generators, and the constraint ensures that the sum of the power that can be generated by the generators is at least the power demand at time t, with high probability of 1 - ϵ. If we consider the worst case (setting ϵ = 0), then the problem can be solved as standard MIP (taking the maximum realization for each random variable ˜ dt). However, the obtained solution is too pessimistic, because it is unrealistic to assume that the randomvariable ˜ dt for power demand attains its maximum value at every t. Therefore it is natural to model the unit commitment problem as CCSP where the random variables appear only in the right-hand side of the constraints. We can easily incorporate complex constraints arising in real schedules of generators into this CCSP formulation.

The unit commitment problem is an important problem in the energy industry and for sustainability, because variants of this problem are also used in large buildings and factories. The large buildings and factories are typically equipped with heat exchanging systems and private power generators. For economical and environmental reasons, it is critical to create an efficient schedule for the cooling machines in the heat exchanging system as well as for the private power generators to satisfy the cold water and power demands of the buildings and factories.

3 Chance-Constrained Stochastic Programming

In this paper, we consider chance-constrained stochastic programming where (i) only the right-hand side vector ˜b is random and (ii) the random vector ˜b has a discrete distribution having a finite support, which can be written as

where c is a vector in Rm, A is an d × m matrix, ˜b = ( ˜ b1, . . . , ˜ bd) is a random vector of values in Rd , ϵ is a parameter between 0 and 1 (typically close to 0), and X is a subset of Rm1 × Zm2 where m1 +m2 = m. The random variables in ˜b may depend on each other. Constraint (2) is a joint probabilistic constraint that ensures that the inequality Ax ≥ ˜b holds with at least probability 1- ϵ. The parameter ϵ determines the tradeoff between the cost and the robustness. That is, a larger ϵ (or equivalently, more robustness) yields a smaller objective function value. We assume that the constraint x ∈ X can be written in the form of the constraints of an MIP and that X is defined by the combination of m1 continuous variables and m2 integer variables.

We assume that each random variable in the vector ˜b = ( ˜ b1, . . . , ˜ bd) has a finite number of realizations. Specifically, the random variable ˜ bi for i ∈ [1, d] takes a realization from hi1, hi2, . . ., hiKi with probability πi1, πi2, . . . πiKi , where ∑k πik = 1. In this paper, for simplicity, we assume that Ki = K for all i. Without loss of generality, we assume that every realization of the random variable ˜ bi takes a non-negative value. We refer to a combination of realizations of d random variables as a scenario, and let S be the set of all scenarios. Specifically, we assume that there are |S| = Kd scenarios h1 , . . . , h|S| , where hs = (hs1 , hs2 , . . . , hsd ) for s ∈ [1, |S|], with probabilities π1, π2, . . . , π|S|, where ∑s∈S πs = 1.

4 SolutionMethods

In Section 4.1, we review a straightforward method to exactly solve the CCSP using MIP, and we improve this method by using a scenario reduction technique in Section 4.2. In Section 4.3, we show another compact MIP formulation to approximately solve the CCSP problem.

4.1 The StandardMethod

Generating all of the scenarios in S explicitly, we can rewrite the CCSP problem in the form of an MIP problem without any random variables.

The binary variable zs is 1 if scenario s does not satisfy Ax ≥ hs . Thus Constraint (7) ensures that the total probability of the relation Ax ≥ ˜b being violated is bounded by ϵ from above. Conversely, zs = 0 implies that under scenario s we have no violated inequality in Ax ≥ hs .

A drawback of this method is the large size of the MIP formulation. Since the number of scenarios |S| is equal to Kd , even if K is small (e.g., K = 3), the number of scenarios is huge (e.g., |S| = 3d ) and exponential in d. Since the size of the MIP grows with increases in |S|, it is unrealistic to exactly solve such a huge MIP.

One of the approaches to tackle this problem is to solve the CCSP approximately by using sample-based techniques, also known as Monte Carlo methods [Campi and Garatti, 2011; Pagnoncelli et al., 2009]. These methods sample a small subset S ′ from S, and solve the MIP with S replaced by S ′ [Luedtke, 2010]. It is known that we require roughly O(|S|/ϵ) samples to yield a highly reliable solution that is feasible and optimal [Luedtke and Ahmed, 2008]. However, it may be computationally prohibitive to solve large problems or to solve problems with high feasibility requirements, because the subset S ′ may be still too large to solve the MIP.

4.2 Reducing the Number of Scenarios

We propose a more efficient approach to exactly solve the CCSP by reducing the number of scenarios, using vector dominance relations. Before showing the details of our approach, we define some notation. Given two vectors ha = (ha1 , . . . , had ) and hb = (hb1 , . . . , hbd ) of length d, the relation ha ≤ hb means that hai ≤ hbi holds for every i ∈ [1, d], whereas the relation ha ̸≤ hb means that there exists at least one i such that hai > hbi .

Definition 4.1 We say that scenario a ∈ S dominates scenario b ∈ S if hb ≤ ha holds. (See Figure 1.)

Definition 4.2 We say that a solution x of a CCSP with probabilistic constraint Pr(Ax ≥ ˜b) ≥ 1 - ϵ is feasible for scenario s if Ax ≥ hs holds.

Definition 4.3 Given a set of scenarios {1, . . . , k} ⊆ S, let h be the vector that is obtained from the elementwise max operation for h1 , . . . , hk . Since the relations hi ≤ h holds for any i ∈ [1, k], we call such a vector h the lowest common ancestor.

Now we have Lemma 4.4.

Lemma 4.4 For any feasible solution x of a CCSP, there exists a scenario s such that (i) x is feasible for any scenario s′ that satisfies hs′ ≤ hs and (ii) x is infeasible for any scenario s′ that satisfies hs′ ̸≤ hs .

Proof. Suppose that the solution x is feasible for a CCSP whose probabilistic constraint is Ax ≥ ˜b. Let Sx be the set of scenarios that are feasible for x. That is, relation Ax ≥ h holds for any scenario h ∈ Sx and relation Ax ̸≥ h holds for any scenario h ̸∈ Sx.

We show that the lowest common ancestor hs of the scenarios in Sx satisfies the conditions (i) and (ii).

Proof for condition (i): We have Ax ≥ hs from the definition of hs . Hence, for any hs′ ≤ hs , we have Ax ≥ hs ≥ hs′ .

Proof for condition (ii): Assume that x is feasible for a scenario s′ such that hs′ ̸≤ hs . Scenario s′ must be in Sx and so hs′ ≤ hs by the definition of hs , which contradicts the assumption. Therefore (ii) holds.

From Lemma 4.4, we know that there exists a scenario s∗ that corresponds to an optimal solution x∗ of a given CCSP instance. Let S ∗ be the set of scenarios that are dominated by s∗ . Then the sum of πs for s ∈ S∗ must be at least 1 - ϵ. Therefore, the set of scenarios that are candidates for the scenario s∗ is given by

By using the set Sϵ, we rewrite the CCSP as

where zs for s ∈ Sϵ is a binary variable. The binary variable zs is 1 if scenario s corresponds to an optimal solution x∗ of the original CCSP. Constraint (14) ensures that one of the scenarios must correspond to the optimal solution x∗ . Conversely, zs = 0 implies that under scenario s we may violate the inequality in the probabilistic constraint (Ax ̸≥ hs ).

Although the size of Sϵ may be close to that of S in the worst case, in practice the size of Sϵ will be far smaller than that of the S. Hence the MIP formulation obtained from this approach is far smaller than the formulation given in Section 4.1, and the solution time of the MIP is significantly reduced.

From Constraint (14), we can see that Constraint (13) can be replaced with

Since Constraint (18) is stronger, in practice we should use this instead of Constraint (13).

4.3 Approximate Solution

As we mentioned in Section 4.2, the number of scenarios Sϵ can still be huge (which means exponential in the number of random variables d, on the order of Kd ). Therefore we consider another approach that approximately solves a CCSP problem. Our new idea is to construct another set S ′ϵ ⊆ Sϵ that can be expressed with a polynomial number of constraints.

Before showing how to construct such a set S ′ϵ , let us define some notation. We say that a random variable ˜ bi touches an undesirable event for scenario s if a realization hik of ˜ bi satisfies the relation hsi ≤ hik. Let Es i be the event in which the random variable ˜ bi touches the undesirable event for scenario s. Then, by definition of Sϵ, scenario s ∈ Sϵ if and only if

Now let us review Boole's inequality (for example, see page 194 in [Comtet, 1974]).

Theorem 4.5

holds for any scenario s. Note that Boole's inequality holds even if the Es i s are dependent. This relation holds with equality when Pr(Es i ∩ Es j ) = 0 holds for every distinct pair of i and j.

Now we define a set of scenarios S ′ϵ such that scenario s ∈ S ′ϵ if and only if

By using Boole's inequality, it is easy to verify that S ′ϵ ⊆ Sϵ. To formulate an MIP with the scenario set S ′ϵ , we transform the following inequality

where ϵi is a constant, into the one having the form suitable for an MIP. To do so, for a fixed i ∈ [1, d], we introduce K binary variables zik ∈ {0, 1} for k ∈ [1,K]. Then an equivalent form of Inequality (19) is given as

because zi1 = · · · = zip = 1 and zi(p+1) = · · · = ziK = 0 holds if and only if scenario s takes a realization such that hsi = hip.

Now we show our approximate MIP model. Whereas the original CCSP problem can be interpreted as finding a scenario s∗ ∈ Sϵ that corresponds to an optimum solution x∗ of the CCSP, we formulate the approximate MIP problem as finding a scenario s ∈ S ′ϵ whose corresponding feasible solution x minimizes the objective value cx. This relation can be modeled as this MIP:

where yi is the i-th element of vector y and ϵi is a continuous variable. Recalling that the combination of Constraints (24) and (25) is equivalent to Constraint (19), it is easy to verify that this MIP model is created as desired. Note that, since Constraint (25) is redundant, we can remove it when we implement this model in practice.

Our approach can be interpreted as dividing the original CCSP with joint probabilistic constraints into a combination of the individual probabilistic constraints. More specifically, Constraint (24) can be interpreted as bounding by ϵi from above the probability that random variable ˜ bi in ˜b violates the relation Ax ≥ ˜b. Therefore Constraint (23) can be interpreted as a combination of d individual probabilistic constraints, and their total violation probability is bounded by ϵ from above. In other words, we can interpret the joint probabilistic Constraint (7) as divided into a set of individual probabilistic Constraints (24), each of which is parameterized by ϵi, using Boole's inequality.

5 Discussion

In this section, we discuss some existing techniques applicable to our approaches.

5.1 Mixing Set Approach

Here we show that we can accelerate solving the approximate CCSP by using the cutting plane method. For the cutting plane method, we exploit an existing technique on the mixing set to generate valid inequalities.

Constraint (22) is the mixing set arising in our approximate CCSP:

where we omit subscript i for simplicity. A stronger (but equivalent) form of this mixing set can be written as

for v such that ∑vk=1 πk ≤ 1 - ϵ and ∑v+1 k=1 πk > 1 - ϵ. G¨unl¨uk and Pochet [2001] studied this mixing set, and obtained many valid inequalities for this mixing set as in Theorem 5.1.

Theorem 5.1 For R = {r1, r2, . . . , ra} ⊆ {1, . . . , v}, the inequalities

where r1 < r2 < · · · < ra and hra+1 = hv+1, are valid for Q and are facet-defining for conv(Q) when r1 = 1.

We can use this result to generate valid inequalities when we solve the CCSP by using an MIP solver (such as IBM ILOG CPLEX) with the cutting plane method. Since the number of inequalities generatable from this theorem are exponential in vi, we cannot generate all of the inequalities in advance, but Luedtke et al. [2010] showed how to generate appropriate inequalities to accelerate solving the MIP in the cutting plane method.

5.2 Logarithmic Number of Binary Variables Approach

Here we reformulate the MIP formulation in Section 4.3 so that the number of binary variables in the formulation is reduced, using an approach similar to the formulation for a piecewise-linear function with a logarithmic number of binary variables [Ibaraki, 1976; Vielma et al., 2010].

Recall that Constraint (22) can be written as the mixing set Q and the inequalities z1 ≥ z2 ≥ · · · ≥ zv always hold in Q. Then we can express Q by using these inequalities

Notice that each zk is not a binary variable but a continuous variable, and that the number of the binary variables is the logarithm of v (≤ K). In these inequalities, each variable αl takes a binary value. Hence only one λk can be 1, and the other λks must be 0. It thus turns out that each zk has a binary value, and the relation z1 ≥ z2 ≥ · · · ≥ zv always holds.

Note that, if λ1 = 1, then all of the zks are 0. If λv = 1, then zv = 0. Recall that, if zv = 1, then all of the zks are 1 and so the chance-constraint of the corresponding random variable cannot be satisfied with probability 1, which yields an infeasible solution as long as ϵ < 1.

The technique presented here is to give theoretical improvement of reducing the number of binary variables. Practical performance would depend on various factors such as the distributions of the random variables in CCSP.

6 Computational Experiments

In this section, we show some experimental results with our new MIP formulations. We compared three approaches:

ALL: theMIP approach that generates all of the scenarios in S as described in Section 4.1,

SMALL: the MIP approach that generates all of the scenarios in Sϵ as described in Section 4.2, where we use the tighter inequality (Constraint (18)), and

APPROX: the compact MIP approach with approximation as described in Section 4.3 using the cutting-plane method mentioned in Section 5.1.

All of the experiments were done on a personal computer, having an Intel Core i5 (M560) running at 2.67 GHz, with 3 GB of RAM using Microsoft Windows XP Professional Service Pack 3. All of the programs were written in Java and compiled with JDK 1.6. We used IBM ILOG CPLEX 12.3 as the MIP solver. All of the execution times are expressed in seconds.

To evaluate our new MIP formulations, we considered a transportation problem similar to [Luedtke et al., 2010]. In this problem, there are a set of suppliers I and a set of customers D. Each supplier i ∈ I has a constant capacity Mi and each customer j ∈ D has a demand, which is a random variable ˜ dj. There is a transportation cost cij (for example, fuel usage for a truck) for shipping each unit of product from supplier i ∈ I to customer j ∈ D. We want to ship the products such that ∑

and

where xij is the amount of product shipped from supplier i to customer j and ϵ is a constant parameter. The objective is to minimize the total shipping cost, which is expressed as

Table 1 summarizes the instances we generated for the transportation problem. The capacity Mi of each supplier i ∈ I is fixed to 1000. The demand ˜ dj of each customer j ∈ D is a discrete distribution having a finite support generated by sampling from the Gaussian distribution with mean 1000 and variance 10. The column 'Realizations' specifies the number of sampling for each demand. The parameter ϵ is set to 0.05.

Table 2 shows computational results for two approaches SMALL and ALL, where each column shows the computation times and each number in parenthesis indicates the number of scenarios generated. These results show that SMALL is clearly faster than ALL for all instances. This is because, whereas ALL generates a large number of scenarios, SMALL generates a very smaller number of scenarios.

Table 3 shows computational results for two approaches SMALL and APPROX, where the fourth column shows the approximation ratio: the objective value of the APPROX divided by the optimal objective value. The number in the parenthesis in the second column shows the number of the scenarios generated in SMALL, and the number in the parenthesis in the third column shows the number of the variables zik generated in APPROX. These results show that both SMALL and APPROX generate very small number of scenarios and the variables zik, and that APPROX is clearly faster than SMALL. This difference comes from the fact that it takes long time to generate the scenarios Sϵ in SMALL (actually, the computation time for solving MIP is less than 0.1 second for every input). Therefore an efficient algorithm for generating Sϵ is required. Finally we see that APPROX exactly solves the instances we used in these experiments (i.e., the approximation ratio is 1.0). We expect that the approximation ratio would become worse (increase) for the instances having more constraints and variables.

7 RelatedWork

The CCSP in which the right-hand side vector is random with a discrete distribution having a finite support has been studied extensively. Ruszczy´nski [2002] formulated the CCSP as an MIP, and some new cutting-plane methods were proposed to solve CCSP problems exactly under the assumption that all of the scenarios happen with equal probability [Luedtke et al., 2010; K¨uc¸ ¨ukyavuz, 2012]. K¨uc¸ ¨ukyavuz [2012] further gives a polynomial-size CCSP formulation under the same assumption.

For general CCSP problems with a discrete distribution having a finite support, Luedtke [2010] gave an algorithm to solve it, extending their algorithm for the CCSP in which the right-hand side vector is random [Luedtke et al., 2010].

8 Conclusion

We showed new approaches to solve the CCSP problems in which only the right-hand side vector is random with a discrete distribution having a finite support, and we showed that these approaches can solve larger problems than the existing approaches. The unit commitment problem is one of the important applications of the CCSP problem, and it is also important for energy savings and sustainability.

An apparent future work is to extend our approaches to general CCSP problems. Another future work is to investigate the tradeoff between the approximation ratio and the computation time of the compact MIP approach. It would also be interesting to investigate whether generalized variants of Boole's inequality (i.e., Bonferroni's inequalities) can be used instead of Boole's inequality.

Acknowledgments

The authors would like to thank anonymous reviewers for their valuable comments.