{'fraction_non_alphanumeric': 0.03586574578940991, 'fraction_numerical': 0.010299285108445414, 'mean_word_length': 4.702326116996776, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chance Constrained Mixed Integer Program: Bilinear and Linear Formulations, and Benders Decomposition

In this paper, we study chance constrained mixed integer program with consideration of recourse decisions and their incurred cost, developed on a finite discrete scenario set. Through studying a non-traditional bilinear mixed integer formulation, we derive its linear counterparts and show that they could be stronger than existing linear formulations. We also develop a variant of Jensen's inequality that extends the one for stochastic program. To solve this challenging problem, we present a variant of Benders decomposition method in bilinear form, which actually provides an easy-to-use algorithm framework for further improvements, along with a few enhancement strategies based on structural properties or Jensen's inequality. Computational study shows that the presented Benders decomposition method, jointly with appropriate enhancement techniques, outperforms a commercial solver by an order of magnitude on solving chance constrained program or detecting its infeasibility.


Introduction
Chance constrained mathematical program (CCMP) first appeared in [7] in 1958 and the (joint) chance constraint was formally introduced in [22,26]. It is often used to capture randomness and to restrict the associated risk on real system operations [32], and is related to optimization problems with Value-at-Risk constraints [31]. Mathematically, it can be represented as the following. Let ω represents a random scenario in Ω, x be the first stage variables and y(ω) be the second stage recourse variables of scenario ω, we have x ∈ R n1 where G(ω) and H(ω) are random matrices, h(ω) is a random column vector of appropriate dimension, and ε is the risk tolerance level. Constraint (3), which is denoted as the chance constraint, enforces that the inequality G(ω)x + H(ω)y(ω) ≥ h(ω) should be satisfied with a probability greater than or equal to 1 − ε. Function F in the objective function (1) represents cost contribution from the recourse decisions. When the model is static such that cost of recourse decisions is ignored, function F is set to 0 [17,29]. Without loss of generality, we assume that the set defined by the first stage constraints is not empty, i.e., X = {x ∈ R n1 + × Z n2 + : Ax ≥ b} = ∅. Because it can explicitly deal with uncertainties and risk requirements, CCMP has been applied to build decision making tools for many real systems where uncertainty is a critical factor in determining system performance. Such applications include service system design and management [11,23], water quality management [36], optimal vaccination planning [38], energy generation [24], and production system design [15]. However, it is well recognized that a CCMP formulation is very challenging to solve, especially when random ω is with a general continuous distribution. In such case, multivariate integration will be involved in and it is extremely hard to derive a closedform expression to represent (3), which makes the application difficult [28]. On the other hand, CCMP with Ω of a discrete and finite support, i.e., Ω = {ω 1 , . . . , ω K }, has received very much attention in recent years. Such CCMP actually is also constructed in sampling based methods to solve those with general continuous distributions [18,6], which makes it practically applicable. Up to now, several advanced and analytical solution methods have been developed for this type of CCMP problems. For the cases where randomness only appears in right-hand-side and there is no recourse decision, the concept of p-efficient point of the distribution [27] is introduced. It leads to a few methods that (partially) enumerate those points to derive optimal solutions [8,3,14]. Because the number of efficient points could be huge, such solution strategy might not be effective in practice. Alternatively, by using a binary indicator variable and "big-M" coefficient, CCMP can be formulated with a set of linear constraints as the following where M is a sufficiently large number.
Note that the enforcement of constraints of scenario ω k is controlled by z k . When z k = 0, those constraints must be satisfied. When z k = 1, they can be ignored due to big-M. Consequently, when F(y) is 0 or linear, CC − bigM is a linear mixed integer program (MIP). Currently, as it can handle more general situations, e.g., recourse decisions are included and the coefficient matrix is random, CC − bigM has become a popular formulation to address CCMP problems [29,16,17,37,35,33]. To improve the computational performance on this type of MIP, a few sophisticated cutting plane methods have been developed [29,37,13,17]. For example, inequalities based on precedence constrained knapsack set are derived and included to improve the solution of the MIP formulation [29]. Valid inequalities based on IIS (irreducibly infeasible subsystem) are derived and implemented within a specialized branch-and-cut procedure [37]. Based on existing research on the mixing set [10], strong inequalities generalizing star inequalities [1] are derived to strengthen the instance only with random right-hand-side [16,13]. Such research is further generalized to compute those with a random coefficient matrix in [17], where strong inequalities that aggregate Benders feasibility cuts are developed using the mixing set structure. The aforementioned algorithms and methods, especially MIP based approaches, have significantly improve our solution capability to compute CCMPs. Nevertheless, those algorithms are complicated for practitioners to use. Also, they often depend on some non-trivial assumptions and are not applicable for solving general type instances, e.g., where cost from recourse decisions should be considered or constraints on recourse decisions could make the whole formulation infeasible. Indeed, compared to the most relevant stochastic programs (SP), which has received numerous research efforts and can be solved by various efficient algorithms in real applications [32,5], solution methods to CCMP are rather very limited and less capable. Given the essential connection between SP and CCMP, we believe that it would be an ideal strategy to make use of the exiting large amount of research results on SP to investigate and develop general and effective algorithms for CCMP, which however receives little attention yet. Hence, in this paper, we seek to make progress toward this end by considering solving CCMP with a set of discrete and finite scenarios.
We begin our exposition in Section 2 by presenting a non-traditional bilinear formulation of CCMP, strong linear reformulations, and Jensen's inequality for CCMP. Inspired by using Benders decomposition method to solve SP, the arguably most popular approach for SP, we then in Section 3 provide a bilinear Benders reformulation, a bilinear variant of Benders decomposition method and a few enhancement strategies, all of which can be easily implemented. In Section 4, we present numerical results obtained from unstructured random instances, which show that the presented bilinear Benders decomposition method typically outperforms the state-of-the-art commercial solver by an order of magnitude, on deriving optimal solutions or reporting infeasibility. Finally, we provide concluding remarks in Section 5.
We mention that our paper provides the following contributions and insights, some of which might be counterintuitive. (i) The bilinear formulation of CCMP is informative. Its linearized formulations can be theoretically stronger and more computationally friendly than the popular CC − bigM formulation. So, it is worth more research efforts to investigate the mathematical structure and to study fast computing algorithms. (ii) A variant of Jensen's inequality is derived for CCMP, which generalizes its SP counterpart. (iii) Different from existing understanding that Benders decomposition may not be a good method [37], Benders decomposition, in the presented bilinear form, is very capable to solve CCMP. Indeed, it is interesting to note that Benders decomposition is probably more effective (in terms of iterations) in solving CCMP than in solving SP. (iv) A few non-trivial enhancement strategies that specifically consider the structure of CCMP are developed, which could lead to further performance improvement. (v) The presented Benders decomposition method yields an easy-to-use fast algorithm strategy that basically does not depend on special assumptions and can solve general CCMPs, which have not been addressed in existing literature. Note that it rather provides a framework that can incorporates numerous existing results on Benders decomposition for SP.

A Bilinear Formulation of CCMP
To simplify our exposition, we define a scenario contributing to the satisfaction of the chance constraint as a responsive scenario, and a scenario that can be ignored (i.e., the corresponding constraints in that scenario can be violated) as a non-responsive scenario. In this paper, we assume that the (expected) recourse cost is derived as the weighted sum of costs from responsive scenarios, which is F(y) = MIP formulation, we present a mixed integer bilinear formulation of CCMP as the following.
x ∈ R n1 The validity of CC − MIBP is obvious. Note that by assigning z k to one or zero, the impact of scenario k, including the cost contribution of recourse decisions and feasibility requirements from recourse constraints, will be removed from the whole formulation or imposed explicitly. Although CC − MIBP is a non-traditional nonlinear formulation, it provides a more direct and natural representation of the underlying combinatorial structure, which captures CCMP without information loss and allows us to gain deep insights. We define θ * (ε) = +∞ if it is infeasible and θ * (ε) = −∞ if it is unbounded. Noting that an instance with a larger ε is a relaxation of another one with a smaller ε, the next result follows.
We mention that when ε = 0, it is clear that z k = 0 for all k and CC − MIBP reduces to a linear MIP, which actually is the underlying SP model. According to Proposition 1, the optimal value of SP, i.e., θ * (0), provides an upper bound to θ * (ε) for any ε > 0.
Remarks: We note that it is often assumed in the study of SP that it has a finite optimal value, which could help researchers focus on solving real life problems that are generally feasible and bounded. Indeed, given the fact that SP is a linear program or mixed integer program, which has many advanced preprocessing techniques developed to detect infeasibility or unboundedness, such assumption might not be restrictive. However, the situation of CC − MIBP (or CCMP in general) is more involved. Different ε may cause CC − MIBP infeasible, finitely optimal, or unbounded. For example, when ε = 0, the resulting SP model must deal with constraints in all scenarios, which could conflict to each other, forcing the whole formulation infeasible. When ε increases to a larger value, we can control binary z variables to remove conflicting scenarios so that constraints from remaining scenarios are compatible, rendering the formulation feasible.
Assumption: In this paper, to help us focus on typical problems, we make an assumption that is independent of ε. We think such independence is of a critical value for real system users as the sensitivity analysis, including feasibility and the total cost, with respect to ε reveals a fundamental feature of CCMP on managing risks and uncertainty. Specifically, we assume that for any scenario ω k , there ∃ x 0 ∈ X such that the recourse problem is finitely optimal, i.e., If such assumption is violated in ω k , there exist three cases which can be processed by the following operations.
In this case, we fix z k = 1 in CC − MIBP. Alternatively, we can update ε = ε − π k and remove constraints and recourse variables of ω k from CC − MIBP.
In this case, we fix z k = 0 in CC − MIBP. Alternatively, we can simply remove constraints and variables of ω k from CC − MIBP. If CC − MIBP with remaining scenarios is feasible, the original formulation is unbounded. Otherwise, both the updated one and original one are infeasible. Case (iii): X = X 1 ∪ X 2 such that η * k (x) = +∞, ∀ x ∈ X 1 and η * k (x) = −∞, ∀ x ∈ X 2 . Basically, such a case is very rare as η k does not have a transition between X 1 and X 2 . Indeed, when X is a feasible set of LP, this situation will not occur. If it does happen, we make use of Branchand-Bound technique by creating two branches for X 1 and X 2 respectively. Then, operations presented for Case (i) and Case (ii) can be applied to those branches for simplification.
Hence, we can intuitively interpret the aforementioned assumption as that it ensures all scenarios behave in a regular fashion because there is no scenario simply causing unboundness or becoming irrevelent. Cases that do not comply with that assumption are rather extreme and they can be dealt with in the preprocess stage. So, our assumption imposes little restriction on CCMP in practice.

Deriving Linear Reformulations from Bilinear Formulation
CC − MIBP's original bilinear form might not be computationally friendly. It would be ideal to convert it into a linear formulation. In particular, we are interested in a tight linear formulation stronger than the popular CC − bigM.
By expanding (8), we have Consider a situation where no recourse decision is involved, i.e., m = 0 in (10) and (7) disappears, and G k = G for all k, i.e., only right-hand-side h k is random. According to [13], we can assume without loss of generality that h k ≥ 0 for k = 1, . . . , K. Note that, no matter which scenarios are responsive in determining an optimal solution, we have Gx ≥ min l=1,...,K {h l }, where min is applied in a component-wise fashion. Hence, projecting out Gxz k through replacing it by min l=1,...,K {h l }z k in (12), we have Because Gx ≥ min k {h k } ≥ 0, it simply leads to the next result.
We can further extend results in Proposition 2 to more general cases where CCMP has recourse decisions. Assume that the set defined by (12) can be equivalently represented by the following linear inequality: where q k is a coefficient vector with an appropriate dimension. Obviously, to ensure the equivalence between (12) (or (8)) and (16), we should have for any x ∈ X, there exists y k such that the inequality G k x + H k y k ≥ q k is satisfied when z k = 1 (noting that the cost contribution of y k is not a concern when z k = 1). Otherwise, additional restriction will be imposed on X, which suggests that the impact of scenario ω k is not removed when z k = 1. Hence, q k,i , the i th component of q k , should have for every i, where Y k represents the set defined by variable bounds of y k . Certainly, setting q k,i to achieve the equality will give us the tightest linear inequality to replace (8). Nevertheless, it involves computing a bilevel robust optimization problem, which is computationally expensive. A simpler strategy is to fix y k to y 0 k ∈ Y k and let q k,i = min{(G k x) i : x ∈ X} + {H k y 0 k } i . Actually, for a situation where 0 ∈ Y k and f k ≥ 0, a strong model only with linear constraints could be derived to represent CCMP.
the following linear model is a valid formulation for CCMP: Proof. For any x, with f k ≥ 0, i.e., recourse decisions incur non-negative cost, it would be optimal to set y k = 0 assuming it is feasible. Under such a situation, it is valid to simplify (7) to η k = f k y k , which ensures there is no cost contribution when z k = 1.
Remarks: (i) Clearly, applying the aforementioned results will produce a strong linear formulation of CCMP, which alleviate our concern on big-M coefficients. Note that deriving q * k,i requires computing an MIP for every k and i, which may incur non-trivial computational expense. One strategy is to relax those MIPs as LPs and derive weaker coefficients with less computational expense. To compute large-scale practical instances with many scenarios and constraints, we believe that computationally more friendly methods are necessary to achieve a balance between the coefficient quality and computational time. In addition, it is worth mentioning that those results do not depend on the variable types of y k for any k. Hence, it can be used to derive strong linear formulations for those whose recourse problems could be mixed integer programs. (ii) Our discussions also show that for a case where q * k,i is not finite for some k and i or cost of some recourse decision is not non-negative, it will be more complicated. We might not be able to derive a finite coefficient for z k in (18), which also implies big-M method might not work under those situations. Hence, a deeper study on those cases is needed to build strong formulations. On this line, an analytical study on deriving strong formulations and valid inequalities is presented in [40].
Instead of using computational methods to derive linear formulations, a simpler and more general approach is to directly linearize CC − MIBP and obtain an MIP in a higher dimension. Given the fact that z k is binary, constraints in (7) and (12) (equivalently (8)) can be easily converted into linear constraints by applying McCormick linearization method on bilinear terms xz k and y k z k [21]. For example, say x j is one component of x. Its bilinear term x j z k (= λ jk ) can be linearized by including the following constraints where U xj is an upper bound of single variable x j . We recognize that such linearization method depends on an upper bound of each individual variable. Nevertheless, it is typically much easier to estimate a strong upper bound for an individual variable, which indicates the linear formulation obtained using McCormick linearization method will be tighter than CC − bigM. Indeed, in many practices, upper bounds of individual variables (e.g., binary or bounded variables) are naturally defined or available [11,37]. In such a case, there is no need to introduce a very big number as the variable upper bound when we perform linearization. Furthermore, this linearization strategy can deal with situations where recourse decisions are with unrestricted cost coefficients, e.g., recourse decisions are with negative cost. So, because of its generality and simplicity, in the remainder of this paper, unless explicitly mentioned, CC − MIBP and its linearized counterpart obtained through McCormick linearization method are interchangeably used. Given the rich and original structural information implied in CC − MIBP, we believe that more mathematical analysis should be done on this bilinear formulation and specialized algorithms need to be developed. Next, following this line, we present a generalization of Jensen's inequality in CCMP and a variant of Benders decomposition customized according to this bilinear form.

Jensen's inequality for CCMP
The combinatorial structure explicitly presented in (7-8) helps us to extend the classical Jensen's inequality [12,19], a strong inequality for stochastic programs, to CCMP. Consider the situation where H k = H and f k = f for all k.
Theorem 4. (Extended Jensen's inequality for CCMP) Let z 0 be an optimal z with respect to x 0 ∈ X and define Eη( , the optimal expected recourse cost from responsive scenarios. If Eη(x 0 ) < +∞, we have Proof. We defineπ k = Furthermore, we defineη(ω) = min{fy : Because H and f are independent of scenarios, we have thatη(ω) is a convex function of random ω. Therefore, following the definition ofπ k , we have kπ kηk ≥ min f y : Also, noting that z 0 is a particular z satisfying K k=1 π k z k ≤ ε, we have Therefore, the desired inequality follows.

Remarks:
When ε = 0, we have z k = 0 for all k in both CC − MIBP and (20). Note that under such situation, CC-MIBP is SP formulation and (20) reduces to which is Jensen's inequality for SP. Nevertheless, as it involves solving an MIP with fractional inequalities, (20) in general is much more complicated than its SP counterpart. Those fractional inequalities can be converted into equivalent bilinear inequalities as follows Hence, with fractional inequalities replaced by (21), the minimization problem in the right-handside of (20) becomes an MIP with bilinear inequalities. Actually, for some special case, those fractional or bilinear inequalities reduce to linear ones.
Corollary 5. Assume π k = 1 K for all k and let L be the integer such that L K ≤ ε < L+1 K . The extended Jensen's inequality in (20) can be simplified as Proof. When π k = 1 K for all k, it is easy to see that there is an optimal z 0 with k z 0 k = L for given x 0 . Hence, based on the proof of Theorem 4, we replace π k by 1 K and replace k z 0 k by L in (20). Then, the inequality in (22) follows.
We mention that Jensen's inequality is not only of theoretical contribution, but also could bring significantly computational benefits. Through replacing x 0 by variable x, the inequality in (20) or (22) can be converted into a valid inequality to bound the optimal expected recourse cost from below. Application of such strategy in SP formulations has been proven to be very successful in computing complicated real-sized problems [9,2]. Hence, it would be interesting to investigate its performance in various CCMP problems.
problem in scenario ω k as follows.
According to our assumption made in Section 2.1 and the duality theory, it is clear that the feasible set of this dual problem is non-empty, regardless of x 0 . Denote all the extreme points of this set by E k = {µ 1 , . . . , µ p k }, all extreme rays by F k = {υ 1 , . . . , υ t k }, and let E and F be the corresponding collections over k. Next, we present Benders reformulation of CC − MIBP.
Theorem 6. (Bilinear Benders Reformulation) The CC − MIBP formulation can be equivalently reformulated as: (i) If all recourse problems are feasible for any x ∈ X, only those in (28), i.e., those defined by extreme points, are needed in BD − MIBP; (ii) if recourse decisions are with zero cost, only constraints in (29), i.e., those defined by extreme rays, are needed in BD − MIBP.
Again, bilinear terms in constraints in (28) and (29) can be linearized using McCormick linearization method, which converts (28) and (29) into linear inequalities and the whole formulation into an MIP. We mention that, comparing to the original CC − MIBP formulation, this Benders reformulation (and its linearized counterpart) has a couple of advantages.

Remarks:
(i) Note that only the first stage decision variables x are involved in this bilinear Benders reformulation and there is no recourse variable. Such an observation actually gives us an essential advantage in eliminating the negative impact of big-M parameters in linearization. As x is typically introduced to model decisions with concrete definitions in practice, e.g., quantity of available resources or binary decisions, their upper bounds can often be obtained from specific applications. Actually, we note that it is straightforward to obtain physically valid upper bounds for the first stage decision variables in the majority of existing CCMP formulations adopted for practical problems. See [24,11,23,34] for a few examples. Therefore, there is no "big-M" issue and our BD − MIBP formulation is strong, which might lead to computational improvements in solving real problems.
(ii) Note also from (28) and (29) (and their linearized counterparts) that z k acts on all constraints generated from extreme points and rays of subproblem/scenario k. Its binary property indicates the inclusion or omission of a whole group of constraints. If z variables are fixed, the resulting problem has a structure same as the Benders reformulation of a stochastic program. Hence, it suggests that, if Benders decomposition method works well for the underlying SP formulation, it should work well for CCMP.
The next result follows directly from Theorem 6 and it provides the basis to design our Benders decomposition method.
Proposition 7. LetÊ k andF k be subsets of extreme points and extreme rays of DP k for all k, andÊ andF be the corresponding collections over k. Then, a partial Benders reformulation with respect toÊ andF , denoted by BD − MIBP(Ê,F ), is a relaxation of BD − MIBP (as well as CC − MIBP), and its optimal value yields a lower bound.
Note that any feasible solution to CC − MIBP provides an upper bound. So, iteratively including constraints from extreme points and rays, i.e., those in (28) and (29) (which are often referred to as Benders optimality and feasibility cuts respectively), as cutting planes could help us generate better lower and upper bounds, yielding a variant of Benders decomposition method based on our bilinear Benders reformulation.
Note that DP k is finitely optimal or unbounded for all k. Hence, if CCMP is unbounded, it must be resulted from min{cx : x ∈ X} = −∞. To avoid such situation, which may cause Benders decomposition method enumerate all extreme points and rays before the observation of unboundedness, we assume that min{cx : x ∈ X} is finite. Let CC − MIBP(x 0 , z 0 ) represent the formulation with x and z fixed to x 0 and z 0 respectively, LB and U B be the current lower and upper bounds, and e be the optimality tolerance. Our procedure, to which we call (bilinear) Benders decomposition, is as follows.

Basic (Bilinear) Benders Decomposition Method
Step 1. (i) If it is infeasible, terminate. Otherwise, obtain its optimal value, V l , and an optimal solution (x 0 , z 0 ).
For all k such that z 0 k = 0, (i) compute subproblem DP k defined in (23)(24)(25) to obtain an optimal solution µ 0 if bounded, or derive an extreme ray υ 0 if unbounded.
Step 3. -Stopping Condition: If UB−LB LB ≤ e, terminate with a solution associated with U B. Otherwise, go to Step 2.
The aforemention algorithm is very close to the standard Benders decomposition method applied to solve SP, except Benders cuts are only generated for scenario/subproblem with z 0 k = 0. We denote this basic implementation as BD0. Actually, such a restriction is not necessary. So, as an alternative approach, we generate Benders cuts from all scenarios/subproblems regardless of z 0 k 's value. We denote it as BD1. As BD1 has almost the same simple structure and complexity as BD0, unless otherwise stated, we call both of them basic Benders decomposition method. Indeed, the minor change made in BD1 leads to a significant computational improvement (see results in Section 4). Note that one possible termination scenario of this Benders decomposition method is the detection of infeasibility specified in Step 2.a.
Next, we present a simple numerical study on a few SP and CCMP instances, which helps us appreciate the connection between SP and CCMP, and provides ideas on improving the aforementioned Benders decomposition procedure in solving CCMP.

Observations and Insights on Computing SP and CCMP
We present a rather small-scale numerical study that involves computing SP formulation, CC − bigM formulation, and the linearized CC − MIBP formulation. Our objective is to learn the connec-tion and the difference between SP and CCMP, the impact of the chance constraint, and to develop a basic understanding on the performance of Benders decomposition for these two types of problems. We simply generate one instance using each of six instance generation combinations described in Section 4. In order to highlight the impact of chance constraint and to minimize the computational challenges from other aspects, we intentionally select instances whose SP formulations can be solved by professional solver CPLEX in a short time (e.g., less than 60 seconds). Overall, we have 6 instances for binary and general integer x, respectively. Our computation environment and related parameters are described in Section 4.
Our results are presented in Table 1-2. In those tables, CPX denotes using CPLEX to compute MIP formulation, BD denotes Benders decomposition method, CPX-bigM represents computing CC − bigM by CPLEX, CPX-MIBP represents computing linearized CC − MIBP by CPLEX. If one instance can not be solved by some method due to time limit, which is set to 3,600 seconds, or memory issue, we use "T" or "M" to represent that entry, respectively, and report available gap information. Row "avg.: XMS" presents the average values over instances whose MIP (big-M, if applicable) formulations can be optimally solved by CPX. To have a fair comparison between SP and CCMP under Benders decomposition method, we select BD1 to compute CCMP, noting that both BD1 and Benders decomposition in SP generate Benders cuts from all scenarios. In our computational study, big-M is set to 10 5 , which is also used as variable upper bound (except binary variables) in McCormick linearization operations.  Clearly, Table 1-2 confirm the well-known result that chance constrained models are more challenging than their SP counterparts, which is more prominent when solver CPLEX is adopted as the solution method. Note that a single chance constraint could easily incur hundreds of times more computational expenses. Through a closer study on Table 1-2, we have a few more interesting observations that have not been fully discussed in previous study. solved to the optimality. Indeed, the bilinear formulation typically can be 2 to 4 times faster than the other one. Hence, we believe that the bilinear formulation could be both theoretically stronger and computationally more friendly. It definitely deserves a further study.
2. Benders decomposition method, implemented in the presented bilinear form, actually displays a very strong capability to compute CCMP, which is different from an understanding made in [37] that it is not a good method. Unlike the solver, which is severely affected by the single chance constraint in every CCMP instance, Benders method often can compute optimal solutions with reasonably more time. According to Table 1, Benders method on average produces an optimal solution with 3-4 times more computational expenses when SP is converted into CCMP. Nevertheless, the professional solver could fail to derive optimal solutions even with hundreds of times more computational expenses. Hence, it can be seen that Bender decomposition method is much more robust to the computational challenge caused by the chance constraint.
3. We further notice two more non-trivial points from Table 1-2 that support the effectiveness of Bender decomposition method. One is that the performance benchmark between Benders decomposition method and the professional solver demonstrates an opposite behavior in SP and CCMP. In SP instances, the solver is typically more efficient. In CCMP instances, however, Benders decomposition method generally performs an order of magnitude faster than the solver. Another one is that, in terms of iterations, CCMP might not be more challenging than SP. For some instances, CCMP formulation may just need one more iterations, e.g., instances B1 and B2, compared to their SP counterparts. Nevertheless, for some other instances, CCMP formulation can be computed with noticeably less iterations, e.g., instances I2 and I3. We think that such observation is counterintuitive and has not been reported in any previous study. It actually indicates that Benders decomposition method has a capability to quickly identify their optimal sets of responsive scenarios from their scenario pools and derive corresponding solutions. Hence, from these two observations, we believe that Benders decomposition method is probably more appropriate to compute CCMP than to compute SP, given that it can effectively deal with the combinatorial structure implied in the chance constraint.
Besides using the presented basic Benders decomposition method, we can design new enhancement strategies or incorporate existing ones developed for solving SP to further improve computational performance. Next, we present a few improvement techniques that make use of the (possible) easiness of the underlying SP formulation or CCMP's structural properties.

Enhancement Strategies of Benders Decomposition
Results in Section 3.2 motivate us to make use of the rather easier SP model to efficiently compute, exactly or approximately, CCMP. Following this line, we present a couple of strategies that can be incorporated into our basic Benders decomposition method in solving CCMP. The first enhancement strategy is to ignore the chance constraint and simply to solve the corresponding SP model by Benders decomposition method to obtain a set of Benders cuts. According to Proposition 7, those Benders cuts, as in (28) and (29), can be employed to initializê E andF in Step 1 of the basic Benders decomposition method. We refer to this strategy as SP based initialization.
Actually, we can directly treat CC − bigM formulation of CCMP as an SP problem, where both x and z are first stage variables, and compute it by Benders method. Nevertheless, the challenge from big-M coefficients of z remains in computing the resulting master problem. Different from the popular idea of big-M concept, we can adopt smaller coefficients for z, to which we denote small-M coefficients, to serve the initialization purpose. Noting that big-M coefficients are larger than small-M coefficients, the next result follows naturally. Proposition 8. The popular CC − bigM formulation of CCMP is a relaxation of the same formulation with small-M coefficients, i.e., the small-M MIP formulation. Therefore, any feasible solution to the small-M MIP formulation is feasible to CC − bigM, and therefore feasible to CCMP model.
It can be seen that the small-M MIP formulation is to approximate CCMP in the primal space, which leads to feasible solutions and strong upper bounds. We can also derive dual information by treating it as a regular SP, whose first stage variables are x and z, and computing it by Benders decomposition method. Then, similar to SP based initialization, the resulting Benders cuts can be employed for initialization. We refer to this strategy as small-M based initialization. As the small-M MIP formulation resembles CC − bigM more than SP does, those Benders cuts probably are more effective than those obtained in SP based initialization. Remark: We note that small-M MIP formulation might provide a flexible framework to analyze and compute CCMP. When M = 0, small-M MIP formulation simply reduces to SP model, which defines a core set inside X that is feasible to all scenarios. It is worth mentioning that core set is convex. When M increases, small-M MIP formulation approaches CCMP, which also expands that core set. Finally, when small-M is sufficiently large, it becomes the popular CC − bigM formulation for CCMP and that set will become the feasible set of CCMP, which is typically non-convex. Actually, this discussion of small-M MIP formulation might reveal the transition between SP and CCMP, noting that SP is easy to solve and CC − bigM is highly difficult to compute. Hence, a trade off between the computational expense and the solution quality could be achieved by adjusting the value of small-M coefficients. Furthermore, we can interpret that both SP and small-M MIP formulations are to compute (approximately) an optimal solution of CCMP and the presented bilinear Benders decomposition is to perform verification of the optimality of that solution.
One more improvement strategy is to employ stronger upper and lower bounds within Benders decomposition procedure. Consider a given z 0 . Let CC − MIBP(·, z 0 ) represent CCMP formulation with z = z 0 , and denote its optimal value byṼ u (z 0 ). Note that CC − MIBP(·, z 0 ) is essentially an SP problem that could be easy to compute. Indeed, if z 0 has been evaluated, we can remove it from master problem to reduce the computational time. It can be achieved by adding the following integer cut |z − z 0 | ≥ 1 into master problem. Given the binary property of z, we can easily simplify it. Let set K 1 be collection of indices of z 0 k = 1 and K 0 be the complement set. The integer cut can be formulated as Accordingly, we modify the basic Benders decomposition procedure to generate stronger bounds. Specifically, let BD − MIBP(Ê,F |C) represent master problem BD − MIBP(Ê,F ) subject to additional conditions C. In Step 1, we include the initialization C = ∅. In Step 2.a, we replace BD − MIBP(Ê,F ) by BD − MIBP(Ê,F |C), and modify Step 2.a.(i) as "If it is infeasible, terminate (with a solution associated with U B if it exists). . . . " In Step 2.c, we replace CC − MIBP(x 0 , z 0 ) by CC − MIBP(·, z 0 ). It is clear that those changes yield new bounds that are stronger than those produced in the basic Benders decomposition procedure. We refer to this improvement as Benders decomposition with integer cuts.
We mention that the presented enhancement strategies are designed based on CCMP structure explicitly. Indeed, the bilinear Benders decomposition method presented in Section 3 is rather a general framework, which can easily be modified to incorporate almost all type of enhancement strategies developed for Benders decomposition in solving SP problems, such as Pareto optimal cut [20,25], multicut aggregation [4,39], and maximum feasible subsystem cut generation [30]. Given the exiting enormous amount of research results on Benders decomposition and SP, it is definitely worth investigating the integration or customization of those results within the presented bilinear Benders variant to deal with complicated real problems. In Section 4, a preliminary study on random CCMP instances is presented to appreciate the basic Benders decomposition method and enhancement strategies described in this section.

Computational Experiments
Computing Environment and Methods: Our computation is made through Concert Technology of CPLEX 12.4, a state-of-the-art professional MIP solver on a Dell Optiplex 760 desktop computer (Intel Core 2 Duo CPU, 3.0GHz, 3.25GB of RAM) with Windows XP platform. Benders decomposition algorithms are implemented in C++ and CPLEX in the same environment. Table 3 summarizes our ten different computing methods. Among them, we mention that CPX represents using CPLEX to compute the currently popular CC − bigM formulation with big-M set to 10 5 . BD2 denotes a BD1 implementation that generates Pareto optimal Benders cuts [20,25], whose effectiveness has been observed in many applications to address SP problems. For BD5-BD8, small-M based initialization is performed with small-M set to 1000. In BD7 and BD8, we adopt a strongest cuts only strategy, which just employs the strongest Benders cut with respect to the current x 0 in each scenario to perform initialization and ignores other Benders cuts derived from computing small-M MIP formulation. In BD8, to generate a pool of strong Benders cut for initialization, we further implement Pareto optimal Benders cuts when computing small-M MIP formulation.
In our computational study, the time limit is set to 3,600 seconds. Also, when computing CC − bigM formulation and master problem of Benders decomposition, the optimality gap e is set to 0.005. For subproblem DP k , the optimality gap e is set to 10 −4 . If initialization operations are involved, we set the optimality gap of master problem in the initialization stage to 0.02 and restrict the total time for the initialization stage less than 500 seconds. Test Bed: Our numerical study concentrates on unstructured instances that are generated randomly. Specifically, let K, I 1 , I 2 n, and m denote the number of scenarios, the number of constraints in (6), the number of constraints in (8), the dimension of x, and the dimension of y k , respectively. In our study, we consider two setups of (I 1 , I 2 , n, m), i.e., T 1 = (10, 30, 20, 40) and T 2 = (20, 50, 30, 70), and three different K values, i.e., 250, 500 and 1000. For each of six resulting combinations, we generate five instances where their coefficients are randomly selected from the ranges described in Table 4. Following those specifications, we create one set of instances with x being binary variables and another set with x being general integer variables, to study algorithm performance with different types of x. For x being general integer variables, we set their upper bounds to 500, which is also used in linearization operations in Benders decomposition. Note that it is not necessary to require x to be integer variables in our algorithm implementation. In order to gain insights on benchmarking Benders decomposition method and the professional solver CPLEX, we primarily consider random instances of a moderate difficulty (i.e., the computational time by CPLEX is more than 500 seconds and has less than 15% optimality gap before time limit when under ε = 0.1). Overall, we have 30 × 2 instances in total. Computational Results: Table 5 and 6 summarize the computational results of 10 different computing methods with ε set to 0.1 for all instances. Column "obj." presents the objective function value of the best feasible solution. Column "sec." records the computational time in seconds. If the algorithm is terminated because of time limit, we label it by "T". If it is terminated due to insufficient memory, we label it by "M". Column "g(%)" keeps the relative gap in percentage when the algorithm terminates before reaching optimality. Column "itr." presents the total number of iterations in Benders decomposition. If any of initialization techniques is applied, the total number of iterations will be displayed as the number of iterations in the initialization stage + that number in the regular computation stage. At the bottom of those tables, we have a few rows to show the average performances of different computing methods. Row "# solved (S)" provides the number of instances solved to optimality. Similarly, row "# unsolved (U)" provides the number of instances unsolved. Row "avg. sec.: S" presents the average computational time in seconds over instances solved to optimality. Row "avg. gap: U" presents the average relative gap in percentage over unsolved instances. As there are a few instances that do not have gap information, we do not include them when computing this row. Row "CPX/BD: XBS" displays the average of computational time ratios between CPLEX and BD, over instances solved by both CPLEX and Benders decomposition. We think those rows supply us critical information to evaluate the computational performance of different methods. Remarks: Based on those two tables, we make a few observations: (i) Benders decomposition method, even in its basic forms BD0 and BD1, shows a significantly better performance than that of the professional solver CPLEX. With the most effective enhancement strategies, Benders decomposition optimally solves 95% of all instances while CPLEX fails to derive optimal solutions for 50% of those instances. For those that can be computed by both CPLEX and Benders decomposition, according to "CPX/BD: XBS", the latter one generally is faster by an order of magnitude. Indeed, given that our testing instances are unstructured, we can anticipate that such improvement can be more prominent in structured real applications where Benders decomposition can be highly customized.
(ii) Although BD1 is obtained by making a simple modification on BD0, its performance is clear superior to BD0. Enhancement strategies, especially SP based and small-M based initialization techniques, could further improve algorithm performance at a significant level. For example, when x are binary, the best performed method is BD5, which is with small-M based initialization. Comparing BD1 and BD5, we note that more than 25% computational time is reduced. And on those solved by CPLEX, it perform 22 times faster. When x are general integers, it is interesting to note that although BD5 probably works faster on those solved by CPLEX, BD3, which is with SP based initialization, performs more stably by solving more instances than BD5.
Other enhancement strategies generally are not as effective as BD3 and BD5. One explanation is that those instances can be computed with a small number of Benders iterations, which may not provide a suitable platform to demonstrate the effectiveness of other strategies. Nevertheless, there may exist some applications where other enhancement strategies could be very effective. (iii) Results on all Benders decomposition implementations with SP or small-M based initialization empirically confirm our discussion following Proposition 8. It is often the case that one single bilinear Benders iteration is needed after the initialization stage is completed (with a high quality feasible solution). Clearly, that single iteration just serves the verification purpose to the optimality of that solution. It is more obvious in BD5, where small-M based initialization is adopted, than BD3, in both Tables. Such observation indicates that this approximation-verification scheme is of a practical interest where the optimality tolerance is an input parameter from system operator. By adjusting the value of M and using this approximation-verification scheme, we will be able to produce a solution satisfying the optimality requirement with a tolerable computational expense for large-scale practical problems.
To study algorithm performance under a different risk tolerance level, we recompute all those instances with ε = 0.05. Results are presented in Table 7, where BD5 and BD3 are selected as Benders decomposition method for instances with binary x and with general integer x, respectively. There are two points worth our attention. First, under ε = 0.05, those instances become significantly easier for CPX. Consequently, the benefit of using Benders decomposition over CPLEX is less attractive. One explanation is that with ε = 0.05, scenario selection is of less freedom, which suggests identifying optimal responsive scenarios is not challenging. So, those instances become more suitable for CPLEX. Another point is the strong capability of detecting infeasibility demonstrated by Benders decomposition method. In Table 7, there are two instances with K = 250 and T 2 that are labeled "inf", i.e., they are infeasible. On the one hand, it takes non-trivial computational expense for CPLEX to report infeasibility. On the other hand, BD5 identifies their infeasibility by a couple of iterations with much less computational time. As we mention in Section 2.1, such strong capability is of a critical value for practitioners to perform sensitivity analysis with respect to ε.
To further explore the solution capabilities of our proposed Benders decomposition method, we generate 30 difficult instances such that they are with general integer x and gap are larger than 15% in CPLEX before time limit. We adopt BD3 as the Benders computing method and present results in Table 8. Again, Benders decomposition method shows a strong computing capability compared to CPLEX. Among all 30 instances, BD3 successfully solves 6 of them. Except two remaining instances where memory limit causes BD3 terminate without gap information, BD3 always ensures a solution that has a significantly smaller gap in all other 22 instances.

Conclusion
In this paper, we study chance constrained mixed integer program with a finite discrete scenario set. We first present a non-traditional bilinear formulation and analyze its structure. Its linear counterparts are derived that could be stronger than popular ones or big-M formulations. We also develop a variant of Jensen's inequality that extends from stochastic program to chance constrained program. To solve this type of challenging problem, we provide a bilinear Benders reformulation and present a bilinear variant of Benders decomposition method. Different from existing understanding that Benders decomposition might not be effective, we observe that the presented implementation method performs drastically better than a state-of-the-art commercial solver by an order of magnitude, on deriving optimal solutions or reporting infeasibility. In addition, a few non-trivial enhancement strategies are designed, which can further improve solution capability of Benders decomposition method. Overall, the presented Benders decomposition method basically does not depend on special assumptions and provides the first easy-to-use fast algorithm to compute general chance constrained program that has not be addressed in existing literature. We mention that the presented bilinear formulation for chance constrained program is informative and its structure is worth of further analysis and investigation. One research direction is to develop strong linear formulations and valid inequalities based on this bilinear representation, especially for those whose recourse problem is an MIP. We also point out that due to its simple scheme, the presented Benders decomposition method is able to integrate almost all improvement strategies and ideas developed for classical Benders decomposition and stochastic programs. So, it is rather an algorithmic framework that can be further extended.
Another future direction is to explore the computational improvement using Jensen's inequality for structured instances arising from real applications. As shown in [9,2], it can successfully improve performance of Benders decomposition to solve sophisticated stochastic programming problems in healthcare applications. So, it would be interesting to study its performance in various chance constrained problems.