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?


Introduction
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 Mathematical Problems in Engineering 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 7 , and Prékopa 8 .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 9 and for a new class of analytical approximation referred to as Bernstein approximations , see the works by Nemirovski and Shapiro 10 .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 17 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 17 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 17 and the other one is an application to the problem of air pollution.

Some Preliminaries: Basic Well-Posedness
In previous work see 17 was introduced the following parametric family of mixed probability stochastic programs P μ, λ : where Q t, λ is the optimal value function of the problem in second stage: min q y : Wy t, y ≥ 0, λ H j y ≥ p j , j 1, . . ., d 2.2 and i H j , j 1, . . ., d, are set-valued mappings from R m to R r with closed graph; ii p j , j 1, . . ., d, are predesigned probability levels; iii if P R s , P R r denote the sets of all Borel probability measures on R s and R r , respectively, we assume that Δ and Λ are subsets of P R s and P R r ; iv C is a close subset of R m .
All remaining vectors and matrices have suitable dimensions.

Mathematical Problems in Engineering 3
This model extends the traditional two-stage linear stochastic program by introducing some probabilistic constraints λ H j y ≥ p j , j 1, . . ., d 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 P μ, λ were studied in 17 .
The major difficulty in understanding the structure of P μ, λ rests in a dilemma about the function Q.
On the one hand, Q is the optimal-value function of a nonlinear program with parameters t and λ, and parametric optimization mainly provides local results about the structure of Q but global results are very scarce and require specific assumptions that are often hard to verify.
On the other hand, Q arises as an integrand in P μ, λ .For studying properties of the related integral we require global information about Q.
From this viewpoint, it is not surprising that most of the structural results about twostage 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 Q we formulate the following general assumptions.
Assumption A.1.For any λ ∈ Λ there exists a nonempty set R λ ⊆ R s and a Lebesgue null set N λ ⊆ R s such that the function Q •, λ is real valued and measurable on R λ , and continuous on where supp μ denotes the smallest closed set in R s with μ-measure one.
Assumption A.3.There exists a real-valued, measurable function h on R s , we call bounding function, with the following properties.

Q-Majorization
It holds that |Q t, λ | ≤ h t for all t ∈ R λ and all λ ∈ Λ.

Integrability
It holds that R s h z μ dz < ∞ for all μ ∈ Δ.

Generalized Subadditivity
There exists a κ > 0 such that h t

Local Boundedness
For each t ∈ R s there exists an open neighborhood of t where h is bounded.
The essence of Assumptions A.1-A.3 is the following: since Q •, λ is the optimalvalue 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 Q •, λ is finite on some set R λ and A.2 Guarantees that the arguments z − Ax are in R λ for all relevant z and x.Otherwise, Q z − Ax, λ 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 G •, μ, λ .
Introducing the exceptional set N λ in Assumption A.1 makes sense, since Q •, λ often lacks continuity on lower-dimensional subsets of its domain of finiteness.
Furthermore, Assumption A.3 ensures an integrable upper bound for the functions |Q • − Ax, λ | when x is varying in some neighborhood.Any other set of conditions ensuring this could be placed instead.
Clearly, h reflects the global growth of |Q •, λ | whose quantitative analysis is acknowledged nontrivial for nonlinear problems.

Applications
Motivated by the study of stochastic programming problems coming from planning and operational management in power generation companies, in previous work see 17 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 x in the model represents generation capacity investment decisions, such as changes continuous of maximum generation capacity for thermal plants, the variable z is a random demand and y 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 pollutionrelated 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 SO 2 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 i 1, 2, . . ., I to be operated over a time horizon with subintervals t 1, 2, . . ., T and a set of control methods j 1, 2, . . ., J. The first stage variables x ijt represent the amount of SO 2 generated from source i, to be mitigated through control measures j in period t under the regulated emission allowance, and c jt is the operating cost of control measure j during period t.The second-stage variables are related to the probabilistic excess SO 2 from source i to be mitigated through control measures j in period t under SO 2 generation rate z ξ , and d jt is the operating and penalty cost for excess SO 2 emission during period t.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 SO 2 abatement. min If we denote by z it ξ the random variable of SO 2 generation rate in source i during period t, the constraints of pollution control demand are x ijt y ijt ξ z it ξ , ∀i, t.

3.2
Finally, the function H y t ξ , ζ represents the accumulation of SO 2 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 i, given the extent j control taken in period t, and the random variable ζ associated with climatic conditions and predicts SO 2 concentrations in a specific area under different meteorological conditions, then we add the probabilistic limitations on the second-stage variables: where p t is the probability levels with which the limitations are to be met.

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: Mathematical Problems in Engineering ξ and ζ represents the independent random variables.
i ξ ∈ Ξ is the possible realizations of the random variable ξ supported on Ξ ⊂ R s .
ii E stands for expectation with respect to the random variable ξ and y ξ ∈ R m for each realization ξ. iii 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 ξ 1 , ξ 2 , . . ., ξ N of N 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 4.9 If we denote by we have that ω y k ≤ 0 is a convex constraints and conservative for each k 1, 2, . . ., N, in the sense that if for Denote by 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 u k ξ k − Ax, for all k 1, 2, . . .N, we have for all k 1, 2, . . ., N and r λ inf q y λ Wy | y ∈ Y .4.17 Note that r λ − λ ξ k − Ax is the dual function corresponding to v, and the master problem can be solved using a differentiable descent method if r λ inf{ q W λ y | y ∈ Y } is strictly concave function over the set {λ | r λ > −∞}.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 v u , 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 v ξ k − Ax , but for each k 1, . . ., N, where λ * k is the Lagrange multiplier associated to linear constraint in the optimal solution of subproblem P K .Therefore, our problem now is how we find specifically the value of this multiplier λ * k .

Normal Distribution
In this section we investigate the case when the random vector Proof.The existence of y * k depends only on whether the feasibility set is not an empty set.By the other hand, ω y is a convex function, and then S is a convex set, so and by continuity of ω, if ω y k > 0, there is α ∈ 0, 1 such that ω y k α 0. Denote by θ α q y k α , then θ α q y * k − y k ≥ 0 for all α ∈ 0, 1 because q y * k ≥ q y k .This implies that the function θ α is monotone increasing in 0, 1 , and therefore θ 0 ≤ θ α ≤ θ 1 5.13 and then we have q y k α ≤ q y * k and y k α ∈ arg min q y | Wy u k , y ≥ 0, ω y ≤ 0 , 5.14 where ω y k α 0. Finally, it is enough to assign to y *

Conclusions
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.

15 Case 1 .|
analyze the two possible cases for each k 1, . . ., N. Let y k ∈ arg min q y | Wy u k , y ≥ 0 .5.If ω y k ≤ 0, then by y * k y k and λ * k is the Lagrange multiplier associated to linear constraint Wy u k .Case 2. If ω y k > 0. Using the results of the previous proposition, we have to find the solution to the penalized problem: min q y C k ω 2 y , s.t.Wy u k , parameter C k 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, γ j is chosen by the limited minimization rule or the Armijo rule, y j k ∈ arg min q 2C k ω y Wy u k , y ≥ 0 5.18 and, if we denote by λ j k the Lagrange multiplier associated to linear constraint Wy u k in 5.18 , and let y * k be the accumulation point of the sequence {y j k }, that is, there is a subsequence {y j k } j∈J that converges to y * k , then we define by λ * k the corresponding limit point of subsequence {λ j k } j∈J .
and W ∈ M s×m R are deterministic matrices and the probability level p ∈ 0, 1 .ivζ∈ Θ is the possible realizations of the random variable ζ supported on Θ ⊂ R r .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 10 and then, have an efficiently solvable deterministic optimization program with the feasible set contained in the chance constrained subproblem.Let H : R m × Θ −→ R, defined by j y , j 1, 2, ..., r are convex, the components ζ j , j 1, 2, ..., r, of the random vector ζ are independent of other random variables and the moment generating functionsM j t : E exp tζ j , j 1, 2, ..., r 4.4are finite valued for all t ∈ R and are efficiently computable.Then, we have that the problem:min q T y | Ax Wy ξ, y ≥ 0 , s.t.: inf t>0 ⎡ ⎣ h 0 y r j 1 tΛ j t −1 h j y − t log p ⎤ ⎦ ≤ 0 4.5is a conservative convex approximation of the chance constrained subproblems 4.2 , for each realizations of the random variable ξ ξ ∈ Ξ ⊂ R s , where 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 10 and the work by Rockafellar and Uryasev 22 .However, our aim in this paper is more focused on showing a numerical methodology to tackle this type of models.