IIS branch-and-cut for joint chance-constrained stochastic programs and application to optimal vaccine allocation

Abstract

We present a new method for solving stochastic programs with joint chance constraints with random technology matrices and discretely distributed random data. The problem can be reformulated as a large-scale mixed 0–1 integer program. We derive a new class of optimality cuts called IIS cuts and apply them to our problem. The cuts are based on irreducibly infeasible subsystems (IIS) of an LP defined by requiring that all scenarios be satisfied. We propose a method for improving the upper bound of the problem when no cut can be found. We derive and implement a branch-and-cut algorithm based on IIS cuts, and refer to this algorithm as the IIS branch-and-cut algorithm. We report on computational results with several test instances from optimal vaccine allocation. The computational results are promising as the IIS branch-and-cut algorithm gives better results than a state-of-the-art commercial solver on one class of problems.

Introduction

In stochastic programming, instead of assuming that all parameter values are deterministically known, a subset of the parameters of a mathematical program are given probability distributions. This paper concerns stochastic programs with joint chance constraints, which can be formulated as follows:

$$\text{SP}:\text{Min}{c}^{\top }x$$

$$\text{s}.\text{t}.\mathit{Ax}⩽b$$

$$\mathbb{P}\left\{T(\tilde{\omega })x⩽r(\tilde{\omega })\right\}⩾\alpha$$

$$x⩾0\text{.}$$

In formulation (1), $x\in {\mathbb{R}}^{n_1}$ is the decision variable vector, $c\in {\mathbb{R}}^{n_1}$ is the cost parameter vector, $A\in {\mathbb{R}}^{n\times m_1}$ is the deterministic constraint matrix, and $b\in {\mathbb{R}}^{m_1}$ is the deterministic right-hand side vector. In this formulation, uncertainty appears as the multi-dimensional random variables $\tilde{\omega }$, that gives rise to the random technology matrix $T(\tilde{\omega })\in {\mathbb{R}}^{n\times m_2}$ and the random right-hand side vector $r(\tilde{\omega })\in {\mathbb{R}}^{m_2}$. Individual outcomes (scenarios) of the random variable are represented as realizations ω ∈ Ω of the sample space. The aim of such a formulation is to find a minimum cost strategy while allowing a subset of the constraints to be violated an acceptable amount of time.

Chance constraints are best suited to optimization problems in which satisfying a certain set of constraints is desirable but may not be done almost surely. For example, it may be too expensive to satisfy the constraints almost surely. Thus chance constraints are used to find optimal solutions such that an acceptable level of reliability is achieved. Applications include telecommunications with companies needing to guarantee a given quality-of-service to their customers [13], air-quality management with the requirement that pollution not exceed prescribed limits too often [2], and also inventory control problems where the upper and lower bounds on stock have to be violated with low probability [15].

In general, stochastic programs with joint chance constraints are hard to solve. The main reason for this is that given a general distribution of the parameter values, the feasible space of the problem is nonconvex [7]. Also, in the case where the random parameters have continuous distributions, evaluating the feasibility of a single candidate solution can require an extremely hard integration [27]. Most research into solution methods for SP (1) has focused either on identifying probability distributions for the parameters with the property that the feasible space is convex, or on methods for solving the problem when the parameters have discrete distributions.

In the case when all the randomness appears in the right-hand side of the formulation, results on distributions that allow the chance constraints to be formulated as convex programs are given in [7], [25]. In the case of discrete probability distributions, most solution methods use integer programming (IP) or other discrete programming techniques. When only the right-hand side is random, solution methods based on a nonconvex reformulation using p-efficient points are given in [13], [12], [29]. IP formulations and solution approaches for the discretely distributed case are addressed in [9], [19]. Another branch of research has been methods based on sampling approximations [18], [21], [22]. These methods use the sample average approximation method [30] to obtain an approximation by replacing the actual probability distribution by an empirical distribution corresponding to a random sample. A different sampling approach, proposed by [8], is to obtain the approximation of a chance-constrained optimization problem by sampling a finite number of its constraints.

When the technology matrix T(ω) is stochastic, problem SP is significantly harder to solve than for the case in which all the uncertainty appears in the right-hand side. With the matrix T(ω) allowed to be continuously distributed, Prekopa [26] presents a few parameter distributions which make the problem convex. In the case where the parameters are joint-normally distributed, An and Eheart [2] present a method for finding upper and lower bounds on the optimal objective function of the problem using the extreme cases of correlation between the random variables. However, there are applications for which these types of distribution assumptions are too stringent. Consequently, there has been a lot of interest in formulations with discretely distributed random parameters, which are often created by sampling.

Given problem SP with decision variables that are pure integer, Tayur et al. [32] solve the problem using an algebraic geometry approach and Aringhieri [3] applies a tabu search for finding feasible solutions. For problems with random T(ω) and possibly continuous decision variables, a deterministic equivalent problem can be formulated as a ‘big-M’ mixed 0–1 integer program [20]. An advantage to such a formulation is that it may be solved directly by a commercial mixed-integer programming (MIP) solver. However, it often results in extremely weak linear programming (LP) relaxations that hinder finding the optimal solution. To strengthen the LP relaxation, Ruszczyński [28] derives cutting planes based on precedence constraint knapsack polyhedra and gives a specialized branch-and-cut algorithm for the problem.

This paper focuses on theoretical results that can be used for general solution techniques for SP. The main significance of these results is that they are valid for joint chance-constrained problems with discretely distributed random technology matrices and right-hand side vectors. We introduce a new class of optimality cuts, called irreducibly infeasible subsystem (IIS) cuts, for strengthening the LP relaxations of the deterministic equivalent MIP reformulation of SP. We also present a method for improving the upper bound found by the algorithm for the case when no IIS cut can be identified. We then derive a branch-and-cut method based on the IIS cuts, termed ‘IIS branch-and-cut’ algorithm’, and discuss its implementation. Finally, we apply the IIS branch-and-cut algorithm to randomly generated large-scale instances arising in optimal vaccine allocation for epidemic prevention.

The rest of the paper is organized as follows: In the next section we give some background on the problem and present the MIP reformulation of the problem that can be solved directly. In Section 3 we derive IIS cuts and an upper bound improvement strategy for the IIS branch-and-cut algorithm. We present and discuss an implementation of the IIS branch-and-cut algorithm in Section 4 and give computational results in Section 5. Finally, we finish with a summary and point further research topics in Section 6.