Research Article | Open Access
A Numerical Method for Two-Stage Stochastic Programs under Uncertainty
Motivated by problems coming from planning and operational management in power generation companies, this work extends the traditional two-stage linear stochastic program by adding probabilistic constraints in the second stage. In this work we describe, under special assumptions, how the two-stage stochastic programs with mixed probabilities can be treated computationally. We obtain a convex conservative approximations of the chance constraints defined in second stage of our model and use Monte Carlo simulation techniques for approximating the expectation function in the first stage by the average. This approach raises with another question: how to solve the linear program with the convex conservative approximation (nonlinear constrains) for each scenario?
Optimization problems involving stochastic models occur in almost all areas of science and engineering. Financial planning or unit commitment in power systems are just few examples of areas in which ignoring uncertainty may lead to inferior or simply wrong decisions.
Stochastic programming models are optimization problems where the decision have to be made under uncertainty because some of the parameters are random variables, and they may use probabilistic constraints and/or penalties in the objective function. In practice, the numerical solvability of the problem plays an important role and there is a tradeoff between correct statistical modeling and computability. For earlier reviews on the various aspects in stochastic programming see, for example, [1–5].
Two-stage stochastic programming is useful for problems where an analysis of strategy scenarios is desired and when the right-side coefficients are random. The main idea of this model is the concept of recourse, which defines possibility to take corrective actions after a realization of the random event. A decision is first undertaken before values of random variables are known, and then, after the random events have occurred and their values are known, a second stage decision is made to minimize “penalties” that may appear because of any infeasibility. For a good introduction and deepen in various aspects of these models, you should see the books in [4, 6].
Chance constrained optimization problems were introduced in Miller and Wagner , and Prékopa . An alternative to the scenario approximation (Monte Carlo sampling techniques) is an approximation based on analytical upper bounding of the probability for the randomly perturbed constraint to be violated. The simplest approximation scheme of this type was proposed in  and for a new class of analytical approximation (referred to as Bernstein approximations), see the works by Nemirovski and Shapiro . Another approximation of probabilities constraints, is by using the Boole-Bonferroni inequalities; see, for example, [11, 12].
When the stochastic program includes nonlinear terms or when continuous random variables are explicitly included, a finite-dimensional linear programming deterministic equivalent does not exist. In this case, we must use some nonlinear programming types of procedures, see, for instance, [3, 13–16].
In previous work (see ) was extended the traditional two-stage linear stochastic program by probabilistic constraints imposed in the second stage. In the next section, we present a summary with assumptions under which the mixed-probability stochastic program is structurally well behaved and stable under perturbation of both probabilities measures. Moreover, in  can be find, under general conditions, first qualitative continuity properties for the expectation of the objective function and the constraint set-valued maps. Hence, we deduced quantitative stability results for the optimal value function and the solution set under perturbations of probabilities measures.
In the third section, two possible applications that could have this model were shown, the first one is a summary of the case of planning and operational management in power generation companies presented in  and the other one is an application to the problem of air pollution.
2. Some Preliminaries: Basic Well-Posedness
In previous work (see ) was introduced the following parametric family of mixed probability stochastic programs : where is the optimal value function of the problem in second stage: and(i), , are set-valued mappings from to with closed graph;(ii), , are predesigned probability levels;(iii)if , denote the sets of all Borel probability measures on and , respectively, we assume that and are subsets of and ;(iv) is a close subset of .
All remaining vectors and matrices have suitable dimensions.
This model extends the traditional two-stage linear stochastic program by introducing some probabilistic constraints , in the second stage of the problem. These types of constraints add nonlinearities to the problem and basic arguments to analyze the well-posedness of were studied in .
The major difficulty in understanding the structure of rests in a dilemma about the function .
On the one hand, is the optimal-value function of a nonlinear program with parameters and , and parametric optimization mainly provides local results about the structure of but global results are very scarce and require specific assumptions that are often hard to verify.
On the other hand, arises as an integrand in . For studying properties of the related integral we require global information about .
From this viewpoint, it is not surprising that most of the structural results about two-stage stochastic programs concern the purely linear and the linear mixed integer cases, that is, the widest problem classes where parametric optimization offers broader results about global stability.
To lay a foundation for the structural analysis of we formulate the following general assumptions.
Assumption A.1. For any there exists a nonempty set and a Lebesgue null set such that the function is real valued and measurable on , and continuous on .
Assumption A.2. It holds that where denotes the smallest closed set in with measure one.
Assumption A.3. There exists a real-valued, measurable function on , we call bounding function, with the following properties. (1) -Majorization
It holds that for all and all .(2) Integrability
It holds that .(3) Generalized Subadditivity
There exists a such that for all .(4) Local Boundedness
For each there exists an open neighborhood of where is bounded.
The essence of Assumptions A.1–A.3 is the following: since is the optimal-value function of a minimization problem it well may attain the values if the problem is infeasible and if the problem is unbounded. Indeed, Assumption A.1 makes sure that is finite on some set and A.2 Guarantees that the arguments are in for all relevant and . Otherwise, would attain infinite values with positive probability, immediately preventing finiteness of the integral: The continuity part of Assumption A.1 together with Assumption A.3 provides a framework for applying dominated convergence to show continuity of .
Introducing the exceptional set in Assumption A.1 makes sense, since often lacks continuity on lower-dimensional subsets of its domain of finiteness.
Furthermore, Assumption A.3 ensures an integrable upper bound for the functions when is varying in some neighborhood. Any other set of conditions ensuring this could be placed instead.
Clearly, reflects the global growth of whose quantitative analysis is acknowledged nontrivial for nonlinear problems.
Motivated by the study of stochastic programming problems coming from planning and operational management in power generation companies, in previous work (see ) was presented an example where was consider a power systems of plants to be operated over a time horizon. In the case of planning and operational management in power generation companies, the first stage variable in the model represents generation capacity investment decisions, such as changes (continuous) of maximum generation capacity for thermal plants, the variable is a random demand and is the second-stage operational variable representing the level of production of energy.
The latter is also limited by emission rights for carbondioxide that may concern single plants or consortia of plants. The level of permitted emission is considered random, since emission rights are traded at predesigned markets via auctions, for instance, whose outcomes are uncertain to market participants. This motivates to model limitations on the operational variables resulting from emission rights by probabilistic rather than deterministic constraints.
Works in [18, 19] have made several applications to model the problem of air pollution; in these papers authers combine different techniques, including two-stage stochastic programming. We now present, based on these previous works, a variation of these models which include restrictions on the type “chance constraints” in the second stage of the model, that is, an example where there are two completely independent probability measures and of different nature.
In air quality management systems, there are uncertainties in a variety of pollution-related processes, such as pollutant characteristics, emission rates, and mitigation measures. These uncertainties would affect the efforts in modeling pollutant. On the other hand, because it is economically infeasible and sometimes technically impossible to design processes leading to zero emission, decision makers and authorities seek to control the emissions to levels at which the effects are minimized. The problem is how to minimize the expected systems cost for pollution abatement while satisfying the policy in terms of allowable pollutant-emission levels.
The SO2 generation rates may vary with the type of coal that is used at the power plants, as well as the related combustion conditions, which could be expressed as a random variable. As an illustrative example, consider a power system consisting of plants to be operated over a time horizon with subintervals and a set of control methods . The first stage variables represent the amount of SO2 generated from source , to be mitigated through control measures in period under the regulated emission allowance, and is the operating cost of control measure during period . The second-stage variables are related to the probabilistic excess SO2 from source to be mitigated through control measures in period under SO2 generation rate , and is the operating and penalty cost for excess SO2 emission during period . In general, it is considered that this cost is much greater than the cost of operating the first stage variables.
The objective is to minimize the total of regular and penalty cost for SO2 abatement. If we denote by the random variable of SO2 generation rate in source during period , the constraints of pollution control demand are Finally, the function represents the accumulation of SO2 in a particular area sensible, such as a city that is surrounded by emission sources or power plants and which depends, on the one hand, on the excess amount of emissions from each source , given the extent control taken in period , and the random variable associated with climatic conditions and predicts SO2 concentrations in a specific area under different meteorological conditions, then we add the probabilistic limitations on the second-stage variables: where is the probability levels with which the limitations are to be met.
4. Numerical Method
In order to have some idea about how the two-stage stochastic programs with mixed probabilities can be treated computationally, we will study the following stochastic linear programming problem: where and represents the independent random variables.(i) is the possible realizations of the random variable supported on .(ii) stands for expectation with respect to the random variable and for each realization .(iii), , and are deterministic matrices and the probability level .(iv) is the possible realizations of the random variable supported on .
The fundamental idea is to gives a convex conservative approximation of the chance constrained subproblems (4.2), for this, we will fallow the work by Nemirovski and Shapiro (see ) and then, have an efficiently solvable deterministic optimization program with the feasible set contained in the chance constrained subproblem.
Let , defined by and we assume that the functions , are convex, the components , , of the random vector are independent of other random variables and the moment generating functions are finite valued for all and are efficiently computable.
Then, we have that the problem: is a conservative convex approximation of the chance constrained subproblems (4.2), for each realizations of the random variable (), where Note that this approximation (it is known as the Bernstein Approximation) is an explicit convex program with efficiently computable constraints and as such is efficiently solvable.
Now, we can use the Monte Carlo simulation, that is, suppose that we can generate a sample of replications of the random vector and then, we can approximate the expectation function by the average and consequently, we have the sample average approximation method: where If we denote by we have that is a convex constraints and conservative for each , in the sense that if for it holds that , then or equivalently and we obtain the following deterministic problem, with nonlinear constraints:
Remark 4.1. It was demonstrated in theoretical studies and numerical experiments that Quasi-Monte Carlo techniques could significantly improve the accuracy of the sample average approximation problem, for a general discussion of Quasi-Monte Carlo methods see the works by Niederreiter in [20, 21]. Moreover, the problem (4.5) is not the only way to get to the conservative convex approximation of the chance constrained problems, we also can use the convex approximation obtained by Conditional Value at Risk (see  and the work by Rockafellar and Uryasev ). However, our aim in this paper is more focused on showing a numerical methodology to tackle this type of models.
Denote by and the convex subsets of the feasible set of problem (4.14) that do not depend on the sample generated by the random vector , then the problem (4.14) can be rewritten as and then, we can take advantage of separability. If we denote where , for all , we have for all and Note that is the dual function corresponding to , and the master problem can be solved using a differentiable descent method if is strictly concave function over the set . However, this last assumption is very restrictive, in fact, in our specific case is not satisfied because the objective function is linear, so we would have to study under what conditions the gradient of the value function , can be explicitly calculated.
Since the master problem (M.P.) has linear constraints, this can be solved using Frank-Wolfe method, for which only we need to know the gradient of , but for each , where is the Lagrange multiplier associated to linear constraint in the optimal solution of subproblem (PK). Therefore, our problem now is how we find specifically the value of this multiplier .
5. Normal Distribution
In this section we investigate the case when the random vector supported on , has all its components normally distributed.
Let us suppose that , , then the moment generating function is defined as for each .
Proposition 5.1. The Bernstein Approximation of the chance constrained subproblems is given by
Proof. As we saw before, the Bernstein Approximation is a conservative convex approximation of the chance constraints defined as and substituting the expression given in (5.2) in the above relationship, we obtain and then Let us denote by the auxiliary function, it is easy to see that the stationary point is a global minimum of the function , therefore, of the given equation (5.6), we can conclude, after some calculations, that and finally, substituting , we have
Proposition 5.2. Let and suppose that . Then, there is such that .
Proof. The existence of depends only on whether the feasibility set
is not an empty set. By the other hand, is a convex function, and then is a convex set, so
and by continuity of , if , there is such that .
Denote by , then for all because . This implies that the function is monotone increasing in , and therefore and then we have and where . Finally, it is enough to assign to .
Now, we analyze the two possible cases for each . Let
Case 1. If , then by and is the Lagrange multiplier associated to linear constraint .
Case 2. If . Using the results of the previous proposition, we have to find the solution to the penalized problem: for a penalty parameter sufficiently large. To resolve this problem, we can apply again the iterative method of Frank and Wolfe, where in each iteration, we solves a linear problem and then, we have where is chosen by the limited minimization rule or the Armijo rule, and, if we denote by the Lagrange multiplier associated to linear constraint in (5.18), and let be the accumulation point of the sequence , that is, there is a subsequence that converges to , then we define by the corresponding limit point of subsequence .
In this paper, we present a strategy or methodology to be followed to solve a two-stage stochastic linear programs numerically, when the chance constraints are included in the second stage. It suggests treating the two measures of probabilities involved in the problem differently. Since the major difficulty of the problem is in the second stage, we chose to assume that we had a sample of replications of random vector involved in the expected value function in the objective function and approximate it by the average. For the case of chance constraints defined in second stage, the main idea was to obtain a convex conservative approximation and then get to an efficiently solvable deterministic nonlinear optimization program for each scenario considered in the previous sample. Since the number of replicas or sample size is generally very large, and for each one must solve a nonlinear optimization problem because a method of decomposition of general deterministic problem were proposed, then although the problem looks very computationally unwieldy for the special case when the random vector of probability constraint of the second stage has all its components normally distributed, an explicitly Bernstein approximation function was obtained and we showed how each nonlinear optimization problem can be solved separately.
The author wish to thank the referees for their careful reading and constructive remark. This work has been supported by CONICYT (Chile) under FONDECYT Grant 1090063.
- J. R. Birge and F. Louveaux, Introduction to Stochastic Programming, Springer Series in Operations Research, Springer, New York, NY, USA, 1997.
- J. Dupačová, “Stability and sensitivity analysis for stochastic programming,” Annals of Operations Research, vol. 27, no. 1–4, pp. 115–142, 1990.
- R. Wets, “Stochastic programming: solution techniques and approximation,” in Mathematical Programming: The State of the Art (Bonn, 1982), A. Bachem, M. Grötschel, and B. Korte, Eds., pp. 566–603, Springer, Berli, Germany, 1983.
- R. J.-B. Wets, “Stochastic programming,” in Optimization, G. L. Nemhauser, A. H. G. Rinnooy Kan, and M. J. Todd, Eds., vol. 1 of Handbooks in Operations Research and Management Science, pp. 573–629, North-Holland, Amsterdam, The Netherlands, 1989.
- R. J.-B. Wets, “Challenges in stochastic programming,” Mathematical Programming, vol. 75, no. 2, pp. 115–135, 1996.
- A. Ruszczyński and A. Shapiro, Eds., Handbooks in Operations Research and Management Science, Volume 10: Stochastic Programming, Elsevier, Amsterdam, The Netherlands, 2003.
- L. B. Miller and H. Wagner, “Chance-constrained programming with joint constraints,” Operations Research, vol. 13, pp. 930–945, 1965.
- A. Prékopa, “On probabilistic constrained programming,” in Proceedings of the Princeton Symposium on Mathematical Programming (Princeton Univ., 1967), pp. 113–138, Princeton University Press, Princeton, NJ, USA.
- A. Ben-Tal and A. Nemirovski, “Robust solutions of linear programming problems contaminated with uncertain data,” Mathematical Programming, vol. 88, no. 3, pp. 411–424, 2000.
- A. Nemirovski and A. Shapiro, “Convex approximations of chance constrained programs,” SIAM Journal on Optimization, vol. 17, no. 4, pp. 969–996, 2006.
- A. I. Kibzun and Y. S. Kan, Stochastic Programming Problems with Probability and Quantile Functions, Wiley-Interscience Serie in Systems and Optimization, Wiley, 1996.
- A. Prékopa, “Boole-Bonferroni inequalities and linear programming,” Operations Research, vol. 36, no. 1, pp. 145–162, 1988.
- J. R. Birge and L. Q. Qi, “Subdifferential convergence in stochastic programs,” SIAM Journal on Optimization, vol. 5, no. 2, pp. 436–453, 1995.
- J. Birge and M. Teboulle, “Upper bounds on the expected value of a convex function using gradient and conjugate function information,” Mathematics of Operations Research, vol. 14, no. 4, pp. 745–759, 1989.
- Y. Ermoliev and R. J.-B. Wets, Eds., Numerical Techniques for Stochastic Optimization, vol. 10 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1988.
- S. Uriasiev, “Adaptive stochastic quasigradient methods,” in Numerical Techniques for Stochastic Optimization, pp. 373–384, Springer, Berlin, Germany, 1988.
- P. Bosch, A. Jofré, and R. Schultz, “Two-stage stochastic programs with mixed probabilities,” SIAM Journal on Optimization, vol. 18, no. 3, pp. 778–788, 2007.
- Y. Li, G. H. Huang, A. Veawab, X. Nie, and L. Liu, “Two-stage fuzzy-stochastic robust programming: a hybrid model for regional air quality management,” Journal of the Air & Waste Management Association, vol. 56, no. 8, pp. 1070–1082, 2006.
- K. Saenchai, L. Benedicenti, and G. H. Huang, “A mixed-integer two-stage interval stochastic programming model for regional air quality management,” Environmental Informatics Archives, vol. 5, pp. 168–176, 2007.
- H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, vol. 63 of CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1992.
- H. Niederreiter, “Quasi-Monte Carlo methods for multidimensional numerical integration,” in Numerical Integration, III (Oberwolfach, 1987), H. Bra and G. Hammerlin, Eds., vol. 85 of International Series of Numerical Mathematics, pp. 157–171, Birkhäuser, Basel, Switzerland, 1988.
- R. T. Rockafellar and S. P. Uryasev, “Optimization of conditional value at risk,” Journal of Risk, vol. 2, pp. 21–41, 2000.
Copyright © 2011 Paul Bosch. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.