Introduction to Stochastic Integer Programming

Shabbir Ahmed

1  Introduction
2  An Example
3  Algorithmic Challenges
4  Theory and Algorithmic Progress
5  Concluding Remarks
6  Links

1  Introduction

This document is part of the Stochastic Programming Community Page (sponsored by the The Committee on Stochastic Programming - COSP) and provides a first introduction to the challenging and exciting field of stochastic integer programming (SIP). Familiarity with basic mathematical programming concepts is assumed. For an introduction to general stochastic programming, please visit the Stochastic Programming Introduction webpage. We illustrate some standard SIP modelling principles using a simple example; discuss the challenges in solving SIP problems; and briefly mention some of the progress in SIP theory and algorithms. For detailed discussion of SIP, the reader is referred to the textbook [5], or the survey articles [13], [16], [21] and [23].

2  An Example

The SunDay Icecream Co. is planning the location and capacity of its distribution centers to serve the demands from its retailers in N cities. The locations and capacities of the distribution centers should be such that retailer demands can be satisfied economically. In addition to the location and capacity decisions for the distribution centers, SunDay also assigns retailers to distribution centers. That is, each retailer's demand is to be supplied from one designated distribution center. This assignment incurs a fixed cost regardless of the retailer's demand.
Unfortunately, the demands of the various retailers are quite uncertain, and SunDay only has information on the probability distribution of the retailer demands. How should SunDay decide on the optimal location and capacity of the distribution centers, as well as the optimal assignment of distribution centers to the retailers? The following subsections describe some alternative optimization models for making this decision.
We shall use the following notation in describing our models. Let yi denote the 0-1 variable indicating whether a distribution center is to be located in city i (i=1,...,N), and xi denote the non-negative variable indicating the capacity of the distribution center located in city i. The maximum capacity that can be located in city i is denoted by Ui. The cost parameters ai and bi denote the per-unit capacity cost and the fixed location cost for locating a distribution center in city i, respectively. Let zij denote the 0-1 variable indicating the assignment of retailer j to distribution center i and gij denote the associated fixed assignment cost. Let ~dj be the uncertain demand of the retailer in city j (j=1,...,N). All costs are assumed to be amortized to a weekly basis, and the demands and capacities are in tonnes of ice-cream.

A simple integer recourse model

Since the demand is uncertain, once the location, capacity and assignment decisions are made, SunDay might find itself in the undesirable situation that the total demand of the retailers assigned to a particular distribution center exceeds capacity of that distribution center. Suppose, in this case, contractual obligations require SunDay to pay a penalty on the shortfall. Furthermore, the penalty cost for a shortfall in a distribution center located in city i is qi dollars per ton of ice-cream shortage rounded up to the next ton. Think of this shortage penalty as the cost to outsource the ice-cream shortfall using one-ton capacity ice-cream trucks.
Note that the locations, capacities and assignments of the distribution centers have to be decided prior to observing the actual demand, whereas the shortage penalty is incurred upon observing the actual demand volume in the future. Consequently, as of now, the shortage penalty is an uncertain parameter, whose actual value depends on the uncertain future demand volume. A natural goal may be to decide on the location, capacity and assignments of the distribution centers here-and-now such that the sum of the location-capacity-assignment costs and the expected future shortage penalty is minimized. An optimization formulation for this problem is then:
[Switch to alternative lay-out]



min

N
å
i=1 
 (ai xi + bi yi) + N
å
i=1 
N
å
j=1 
 gijzij + E[Q(x,z,
~
d
 
 )]
s.t.
0 £ xi £ Ui yi           i=1,...,N


N
å
i=1 
 zij = 1            j=1,...,N


N
å
j=1 
 zij £ yi           i=1,...,N

yi,  zij Î {0,1}         i,j = 1,...,N

where x : = (x1,...,xN)T, z : = (z11, z12, ...,zNN)T, ~d : = (~d1,...,~dN)T,
Q(x,z,
~
d
 
 ) : = N
å
i=1 
 qi é
ê
max
 { N
å
j=1 
~
dj
 
 zij - xi, 0} ù
ú 		,
and éa ù denotes the ceiling or integer round-up of a.
The above problem is an example of a two-stage stochastic program with simple integer recourse. The first-stage decisions constitutes location-capacity-assignment of the distribution centers prior to observing the exact demands, while the second-stage decisions constitute a simple recourse action of paying a penalty on the shortfall once demand information becomes available. Since the penalty is paid on integer shortfall values, the recourse action is "integer."

A general integer recourse model

Suppose now that SunDay has the flexibility of postponing the assignment decisions until actual demand information becomes available. In this case, a natural model is to minimize the sum of location-capacity costs and the expected future assignment and shortage penalty costs. The formulation becomes


min

N
å
i=1 
 (ai xi + bi yi) + E[Q(x,
~
d
 
 )]
s.t.
0 £ xi £ Ui yi           i=1,...,N

yi Î {0,1}               i = 1,...,N

where

Q(x,
~
d
 
 ) : =

min

N
å
i=1 
N
å
j=1 
 gijzij + N
å
i=1 
 qi wi

s.t.

N
å
j=1 
~
dj
 
 zij - wi £ xi       i=1,...,N



N
å
i=1 
 zij = 1                  j=1,...,N


wi Î Z+,  zij Î {0,1}     i,j = 1,...,N

The above problem is an example of a two-stage stochastic program with general integer recourse. Here the second-stage or recourse action constitute optimally deciding the assignment of the retailers to distribution centers and computing any shortages based on the capacity installed in the first-stage and actual demand realized.

A chance-constrained model

Suppose that, as in the simple integer recourse example, SunDay has to make the location-capacity-assignment decisions prior to observing demand. However, instead of paying a penalty on shortages, it wishes to guarantee that the probability of a shortage at any of the distribution centers is small ( £ e). The problem formulation now becomes:


min

N
å
i=1 
 (ai xi + bi yi)+ N
å
i=1 
N
å
j=1 
 gijzij
s.t.
0 £ xi £ Uiyi                   i=1,...,N


N
å
i=1 
 zij = 1                   j=1,...,N


N
å
j=1 
 zij £ yi           i=1,...,N

yi,  zij Î {0,1}              i,j = 1,...,N


Pr
    æ
è 		N
å
j=1 
~
dj
 
 zij £ xi,  i=1,...,N ö
ø 		  ³ 1 - e

The last constraint in the above model bounds the probability of a shortage from above. Such constraints are known as chance constraints or probabilistic constraints. Unless the underlying random parameter distributions and/or the constraint structure is of very specific form, chance-constraints are extremely difficult to deal with algorithmically. For a discussion of general chance-constrained stochastic programs, see [18]. Various special classes of chance-constrained SIPs have been studied, see e.g. [3] and [4].
Owing to the difficulties with chance-constrained SIPs, much of the progress in SIP theory and algorithms has been for recourse models. For the remainder of this article, we shall concentrate on recourse SIP models.

[Top of page]

3  Algorithmic Challenges

In this section, we discuss some of the computational challenges involved in solving stochastic programs with integer recourse. A standard form for a two-stage stochastic program with integer recourse is as follows:


min
cTx + E[Q(x,w)]
s.t.
Ax = b,

x Î R+n1-p1 ×Z+p1

where

Q(x,w) : =

min
qTy

s.t.
Wy = h - Tx,


y Î R+n2-p2 ×Z+p2.

Above x represents the first-stage decisions and y represents the second-stage decisions, and w represents the uncertain data for the second-stage (the parameters (q,W,h,T) are actual realization of the random data). Note that both first- and second- stage variables are restricted to be mixed-integer (R and Z denotes reals and integers respectively).
There are three levels of difficulty in solving stochastic integer programs of the above form.
  1. Evaluating the second-stage cost for a fixed first-stage decision and a particular realization of the uncertain parameters. Note that this involves solving an instance of the second-stage problem which may be an NP-hard integer program and involve significant computational difficulties.

  2. Evaluating the expected second-stage cost for a fixed first-stage decision. If the uncertain parameters have continuous distribution, this involves integrating the value function Q(x,·) of an integer program, and is in general impossible. If the uncertain parameters have a discrete distribution, this involves solving a (typically huge) number of similar integer programs.

  3. Optimizing the expected second-stage cost. It is well known that the value function of an integer program is non-convex and often discontinuous. Consequently, the expected second-stage cost function E[Q(·,w)] is non-convex in x. Figure 1 illustrates the non-convex nature of the objective function cTx + E[Q(x,w)] for a small stochastic integer programming problem with two first-stage variables. The optimization of such a complex objective function poses severe difficulties.

newval.png
Figure 1: Objective function of a small SIP

[Top of page]

4  Theory and Algorithmic Progress

In this section, we briefly mention some of the theoretical and algorithmic progress towards addressing the afore-mentioned difficulties in two-stage stochastic integer programming.

Difficulty 1

It is typically assumed that a single evaluation of the second-stage problem is somehow tractable. Without this assumption, very little progress in optimizing the expected value of this integer program is possible. There has been some work in obtaining approximate solutions to SIPs through approximate solutions to the integer second-stage problem, e.g. in [9].

Difficulty 2

Let us now consider the difficulty of evaluating the expected value of the second-stage integer program E[Q(x,w)] for a given first-stage decision x. As mentioned earlier, if the distribution of the uncertain parameters is continuous or if, in case of discrete distributions, the number of possible realizations is extremely large, then it is practically impossible to evaluate E[Q(x,w)] exactly. In this case, one has to resort to approximating the underlying probability distribution by a manageable distribution.
For example, if the underlying distribution is continuous one may approximate it via discretization. Theoretical stability results for SIPs (see [19,20]) suggest that if the approximate distribution is not too "far" from the true distribution, then the optimal solution to the SIP involving the approximate distribution is close to the true optimal solution.
Alternatively, one may use statistical estimates of the expected value function via Monte Carlo sampling. This can be done in one of two ways. In interior sampling approaches, the estimation of E[Q(x,w)] is carried within the algorithm used to optimize this function. For example, in the stochastic branch and bound algorithm [17], the feasible domain of the first-stage variables x is recursively partitioned into subsets, and statistical upper and lower bounds on the objective function cTx + E[Q(x,w)] over these subsets are obtained via sampling. These bounds are used to discard inferior subsets of the feasible domain, and further partition the promising subsets to eventually isolate a subset containing an approximate optimal solution. In exterior sampling approaches, the sampling and optimization are decoupled. A Monte Carlo sample of the uncertain parameters is generated, and the expectation objective in the problem is replaced by a sample average. The resulting approximation of the problem is then solved, and its solution serves as a candidate solution to the true problem. By repeating the sampling-optimization procedure several times, it is possible to obtain statistical confidence intervals on the obtained candidate solutions. Surprisingly, it can be shown that the number of samples needed to get a fairly accurate solution with high probability is quite small. Discussion of theoretical and algorithmic issues pertaining to the above approach in the context of SIP can be found in [14,1].
Regardless of how the underlying distribution is approximated, an evaluation of the expected second-stage objective value (under the approximate distribution) requires solving many similar integer programs. Owing to the absence of a computationally useful duality theory for integer programming, it is very difficult to take advantage of the similarities in the different second-stage IPs. When the second-stage variables are pure integer, several proposals for using Groebner basis and other test set based methods from computational algebra for exploiting IP problem similarity have been put forth [10,22,26]. For the case of mixed-integer subproblems, if a cutting plane method is used, then under some conditions it is possible to transform a cut (or a valid inequality) derived for one of the second-stage subproblems into a cut for another subproblem by exploiting similarity [6,11].

Difficulty  3

Much of the development in SIP has been towards the difficulty of optimizing f(x) : = cTx + E[Q(x,w)], i.e., the sum of the first-stage and the expected second-stage costs. We classify these developments as follows.

Convex approximations of the value function

In case of SIP with simple integer recourse, a single evaluation of f(x) is easy, however owing to the non-convex nature of E[Q(x,w)], the function f(x) is difficult to optimize. Fortunately, it has been shown [12] that an expectation of the continuous counterpart of the simple integer recourse function Q(x,w) (i.e. where the recourse variables are continuous instead of being discrete valued) under a particular class of distributions of w provides the tightest convex under-estimator for E[Q(x,w)] over its entire domain. Recently, similar results for constructing convex approximations of general integer recourse functions (SIPs involving pure integer second-stage variables) by perturbing the underlying distribution have been obtained [27]. These convex approximating functions are amenable for optimization and can be used to provide strong lower bounds within some of the under-mentioned algorithms for optimizing f(x).

Stage-wise decomposition algorithms

This class of algorithms adopt the natural viewpoint of optimizing the objective function f(x) : = cTx + E[Q(x,w)] over the set of feasible first-stage decisions (say denoted by X). Note that E[Q(x,w)] is not available in closed-form, nor is it suited for direct optimization. Typical algorithms in this class progress in the following manner. In an iteration i, the algorithm builds and/or refines a computationally tractable approximation (typically an under-estimator) ÙQi(x) of E[Q(x,w)]. The under-estimating function cTx +ÙQi(x) is optimized with respect to the first-stage variables (this optimization problem is often referred to as the master problem) to obtain a lower bound on the true optimal objective value as well as to provide a candidate first-stage solution xi. Corresponding to the candidate solution, the second-stage expected value function E[Q(xi,w)] is evaluated. Assuming that the distribution of w is discrete, this step involves independent solution of the second-stage problems for each realization of w, allowing for a computationally convenient decomposition. The value cTxi +E[Q(xi,w)] provides an upper bound on the optimal objective value. The evaluation of E[Q(xi,w)] also provides information on how the approximation ÙQi is to be updated/refined to ÙQi+1 for the master problem of iteration i+1. The process continues until the bounds have converged. The details of the various stage-wise decomposition algorithms differ mainly in how the approximation ÙQi is constructed and updated.
For SIPs with binary first-stage variables and mixed-integer second-stage variables, the integer L-shaped method [15] approximates the second-stage value function by linear "cuts" that are exact at the binary solution where the cut is generated and are under-estimates at other binary solutions. The optimization of the master problem, i.e. optimizing cTx + ÙQi(x) with respect to the first-stage binary variables, is carried out via a branch-and-bound strategy. As soon as a new first-stage binary solution is encountered in the branch-and-bound search, the second-stage subproblems are solved to generate a new cut and to refine the approximation ÙQi. The integer L-shaped method requires that the second-stage integer problems (corresponding to the current candidate first-stage solution xi) are all solved to optimality before a valid cut can be generated. Recall that typical integer programming algorithms progress by solving a sequence of intermediate linear programming problems. Using disjunctive programming techniques, it is possible to derive cuts from the solutions to these intermediate LPs that are valid under-estimators of E[Q(x,w)] at all binary first-stage solutions [11,24]. This avoids solving difficult integer second-stage problems to optimality in all iterations of the algorithm, providing significant computational advantage.
For SIPs where the first-stage variables are not necessarily all binary, dual functions from the second-stage integer program can, in principle, be used to construct cuts to build the approximation ÙQi [8]. Owing to the non-convex nature of IP dual functions, the cuts are no longer linear, resulting in a non-convex master problem. If the second-stage variables are pure integer (and the first-stage variables are mixed-integer), then it can be shown that E[Q(x,w)] is piece-wise constant over subsets that form a partitioning of the feasible region of x [22]. Optimization of cTx +E[Q(x,w)] over such a subset is easy. This leads immediately to a scheme where the subsets are enumerated, and the one over which the objective function value is least is chosen. By exploiting certain monotonicity properties, the subsets can be enumerated efficiently within a branch-and-bound strategy [2].

Scenario-wise decomposition

Assuming the distribution of w is discrete, i.e. the random parameter takes one of a finite set of values (scenarios) {w1,...,wS} having probabilities {p1,...,pS}, the two-stage SIP can be re-formulated as follows


min

S
å
s=1 
 ps (cTxs + qsTys)
s.t.
Axs = b   s=1,...,S,

Tsxs + Wsys = hs   s=1,...,S,

xs Î R+n1-p1 ×Z+p1   s=1,...,S,

ys Î R+n2-p2 ×Z+p2   s=1,...,S,

x1 = x2 = ¼ = xS.

Note that copies of the first-stage variable have been introduced for each scenario. The last constraint, known as the non-anticipativity constraints guarantee that the first-stage variables are identical across the different scenarios. Consider the Lagrangian dual problem obtained by relaxing the non-anticipativity constraints through the introduction of Lagrange multipliers. Note that for a given set of multipliers, the problem is separable by scenarios, thus the dual function can be evaluated in a decomposed manner. Optimization of the dual function can be performed using standard non-smooth optimization techniques. However, owing to the non-convexities, there exists a duality gap, and one needs to resort to a branch-and-bound strategy to prove optimality [7].

[Top of page]

5  Concluding Remarks

This article offers a very limited view of some of the general modelling and algorithmic concepts in stochastic integer programming. The concepts alluded to here have been significantly enriched and extended in recent years. One particularly important area not discussed here is the multi-stage extension of two-stage SIP. We have also not mentioned the large number of important developments in application-specific areas of SIP (see, e.g., [25] for a bibliography of applications of SIP). We hope that this simple introduction will pique the readers interest towards further exploration of SIP.

[Top of page]

6  Links

Reference Material

Bibliography

Test Problems

[Top of page]

References

[1]
S. Ahmed and A. Shapiro. The sample average approximation method for stochastic programs with integer recourse. Optimization-online, available at http://www.optimization-online.org , 2002.
[2]
S. Ahmed, M. Tawarmalani, and N. V. Sahinidis. A finite branch and bound algorithm for two-stage stochastic integer programs. Stochastic Programming E-Print Series, available at http://www.speps.info , 2000.
[3]
P. Beraldi and A. Ruszczy\'nski. A branch and bound method for stochastic integer problems under probabilistic constraints. Optimization Methods and Software, 17(3):359-382, 2002.
[4]
P. Beraldi and A. Ruszczy\'nski. The probabilistic set-covering problem. Operations Research, 50(6):956-967, 2002.
[5]
J. R. Birge and F. Louveaux. Introduction to Stochastic Programming. Springer, New York, NY, 1997.
[6]
C. C. Caroe. Decomposition in stochastic integer programming. PhD thesis, University of Copenhagen, 1998.
[7]
C. C. Caroe and R.Schultz. Dual decomposition in stochastic integer programming. Operations Research Letters, 24:37-45, 1999.
[8]
C. C. Caroe and J. Tind. L-shaped decomposition of two-stage stochastic programs with integer recourse. Mathematical Programming, 83:451-464, 1998.
[9]
M. A. H. Dempster, M. L. Fisher, L. Jansen, B. J. Lageweg, J. K. Lenstra, and A. H. G. Rinnooy Kan. Analytical evaluation of hierarchical planning systems. Operations Research, 29:707-716, 1981.
[10]
R. Hemmecke and R. Schultz. Decomposition of test sets in stochastic integer programming. Mathematical Programming, 94:323-341, 2003.
[11]
J. L. Higle and S. Sen. The C3 theorem and a D2 algorithm for large scale stochastic integer programming: Set convexification. Working paper, University of Arizona, 2000.
[12]
W. K. Klein Haneveld, L. Stougie, and M. H. van der Vlerk. On the convex hull of the simple integer recourse objective function. Annals of Operational Research, 56:209-224, 1995.
[13]
W. K. Klein Haneveld and M. H. van der Vlerk. Stochastic integer programming: General models and algorithms. Annals of Operational Research, 85:39-57, 1999.
[14]
A. J. Kleywegt, A. Shapiro, and T. Homem-De-Mello. The sample average approximation method for stochastic discrete optimization. SIAM Journal of Optimization, 12:479-502, 2001.
[15]
G. Laporte and F. V. Louveaux. The integer L-shaped method for stochastic integer programs with complete recourse. Operations Research Letters, 13:133-142, 1993.
[16]
F.V. Louveaux and R. Schultz. Stochastic Integer Programming. In Handbooks in Operations Research and Mangament Science 10: Stochastic Programming. A. Ruszczy\'nski and A. Shapiro (Eds.). North-Holland, 2003 (To appear).
URL: http://www.uni-duisburg.de/FB11/PUBL/SHADOW/558.rdf.html
[17]
V. I. Norkin, G. C. Pflug, and A. Ruszczy\'nski. A branch and bound method for stochastic global optimization. Mathematical Programming, 83:425-450, 1998.
[18]
A. Prekopa. Stochastic Programming. Kluwer Academic Publishers, Netherlands, 1995.
[19]
R. Schultz. Continuity properties of expectation functions in stochastic integer programming. Mathematics of Operations Research, 18(3):578-589, 1993.
[20]
R. Schultz. On structure and stability in stochastic programs with random technology matrix and complete integer recourse. Mathematical Programming, 70(1):73-89, 1995.
[21]
R. Schultz, L. Stougie, and M. H. van der Vlerk. Two-stage stochastic integer programming: A survey. Statistica Neerlandica. Journal of the Netherlands Society for Statistics and Operations Research, 50(3):404-416, 1996.
(Online version)
[22]
R. Schultz, L. Stougie, and M. H. van der Vlerk. Solving stochastic programs with integer recourse by enumeration: A framework using Gröbner basis reductions. Mathematical Programming, 83:229-252, 1998.
[23]
S. Sen. Algorithms for stochastic mixed-integer programming models. Technical report, MORE Institute, SIE Department, University of Arizona, Tucson, AZ, 2003.
URL: http://tucson.sie.arizona.edu/MORE/papers/SIPHbook.pdf
[24]
H. D. Sherali and B. M. P. Fraticelli. A modification of Benders' decomposition algorithm for discrete subproblems: An approach for stochastic programs with integer recourse. Journal of Global Optimization, 22:319-342, 2002.
[25]
L. Stougie and M. H. van der Vlerk. Stochastic integer programming. In Annotated Bibliographies in Combinatorial Optimization, M. Dell'Amico et al (eds), John Wiley & Sons, New York, pages 127-141, 1997.
[26]
S. R. Tayur, R. R. Thomas, and N. R. Natraj. An algebraic geometry algorithm for scheduling in the presence of setups and correlated demands. Mathematical Programming, 69(3):369-401, 1995.
[27]
M.H. van der Vlerk. Convex approximations for complete integer recourse models. Technical Report 02A21, SOM, University of Groningen, The Netherlands, 2002. Also: Stochastic Programming E-Print Series, available at http://www.speps.info
[Top of page]


File translated from TEX by TTH, version 3.49.
On 13 Feb 2004, 18:27.