A smoothing SAA method for a stochastic mathematical program with complementarity constraints

A smoothing sample average approximation (SAA) method based on the log-exponential function is proposed for solving a stochastic mathematical program with complementarity constraints (SMPCC) considered by Birbil et al. (S. I. Birbil, G. Gürkan, O. Listes: Solving stochastic mathematical programs with complementarity constraints using simulation, Math. Oper. Res. 31 (2006), 739–760). It is demonstrated that, under suitable conditions, the optimal solution of the smoothed SAA problem converges almost surely to that of the true problem as the sample size tends to infinity. Moreover, under a strong second-order sufficient condition for SMPCC, the almost sure convergence of Karash-Kuhn-Tucker points of the smoothed SAA problem is established by Robinson’s stability theory. Some preliminary numerical results are reported to show the efficiency of our method.


Introduction
Our concern in this paper is to investigate the almost sure convergence of a numerical method for the following stochastic mathematical program with complementarity constraints (SMPCC): min E[f (x, y, ξ(ω))] (1.1) where Ψ(x, y) := E[F (x, y, ξ(ω))], F : R n × R m × R k → R m is a random mapping, ξ : Ω → Ξ ⊂ R k is a random vector defined on a probability space (Ω, F , P), E denotes the mathematical expectation. Throughout the paper, we assume that E[f (x, y, ξ(ω))] and E[F (x, y, ξ(ω))] are all well defined and finite for any (x, y) ∈ R n × R m . To ease the notation, we write ξ(ω) as ξ and this should be distinguished from ξ being a deterministic vector of Ξ in a context.
Over the past several decades, the mathematical program with complementarity constraints (MPCC) has been intensively studied for its extensive application in engineering, economics, game theory and networks, see [10] and [11] for systematic expositions, examples, and applications. While in practice, there are some important instances where problem data contain some uncertain factors, and consequently various formulations of SMPCC models are proposed to reflect the uncertainties. See [26], [22], [9], [2], [12], and references therein. Among these formulations, Birbil et al. [2] were the first to treat SMPCC (1.1) and presented the so-called sample-path method for solving it.
Evidently, if the integral involved in the mathematical expectation of problem (1.1) exists or is computable, then problem (1.1) is reduced to the usual MPCC problem. However, in many cases, exact evaluation of the expected value in (1.1) is either impossible or prohibitively expensive. Sample average approximation (SAA) method, also known as sample path optimization (SPO) method [18], is suggested by many authors to handle this difficulty, see the recent works [15], [23], [6] and a comprehensive review by Shapiro [24]. The basic idea of SAA is to generate an independent identically distributed (iid) sample ξ 1 , . . . , ξ N of ξ and then approximate the expected value by sample average. In this context, let ξ 1 , . . . , ξ N be an iid sample, then the SMPCC (1.1) is approximated by the following SAA problem: f (x, y, ξ i ) is the sample-average function of f (x, y, ξ) and F (x, y, ξ i ) is the sample-average mapping of F (x, y, ξ). We refer to (1.1) as the true problem and (1.2) as the SAA problem to (1.1).
In this paper, we propose a smoothing SAA method based on the log-exponential function [20] for solving (1.1). That is, utilizing the properties of the log-exponential function, we reformulate the SAA problem (1.2) as a smooth nonlinear programming (NLP) problem by displacing the difficult equilibrium constraints of (1.2) with a smooth function, and then solve the true problem (1.1) by solving a sequence of such smoothed SAA problem.
Recently, the smoothing SAA method has caught great attention among researchers in SMPCC, see e.g., Shapiro and Xu [22], Xu [26], Xu and Meng [27]. However, most available results discuss an application of the SAA method to SMPCC with the assumption that the constraint of (1.1) has a unique solution for every x. Our paper, without such assumption, focuses on the sufficient conditions ensuring the almost sure convergence of the proposed smoothing SAA method. Also notice that although the idea of the SAA method is essentially the same as that of the sample-path method, our method differs from the work of Birbil et at. [2] which discusses the almost sure convergence of the optimal solutions of SAA problem (1.2) without referring to a particular smoothing function.
By the notion of epi-convergence in [20], we establish the almost sure convergence of the optimal solution of a smoothed SAA problem as the sample size tends to infinity. As it is more practical to find a stationary point than a global minimizer of the true problem, under the MPCC strong second order sufficient condition (MPCC-SSOSC) in [21], we investigate the almost sure convergence of the stationary points of smoothed SAA problem.
This paper is organized as follows: Section 2 gives preliminaries needed throughout the whole paper. In Section 3, by introducing an iid sample and the log-exponential function, we formulate the SAA problem (1.2) as a smooth NLP problem. In what follows, we discuss the almost sure convergence of optimal solutions and stationary points of the smoothed SAA problem as the sample size tends to infinity in Section 4. In Section 5, we report some preliminary numerical results.

Preliminaries
Throughout this paper we use the following notation. Let · denote the Euclidean norm of a vector or the Frobenius norm of a matrix. For an m × n matrix A, A ij denotes the element of the ith row and jth column of A. We use B to denote the closed unit ball and B(x, δ) the closed ball around x of radius δ > 0. For an extended real-valued function ϕ : R n → R ∪ {±∞}, epi ϕ, ∇ϕ(x), and ∇ 2 ϕ(x) denote its epigraph i.e. the set {(x, α) : ϕ(x) α}, the gradient of ϕ at x, and the Hessian matrix of ϕ at x, respectively. If a mapping F : R n × R m → R m is Fréchetdifferentiable at (x, y) ∈ R n × R m , then we use J F (x, y) and J x F (x, y) to denote respectively the Fréchet-derivative of F at (x, y) and the partial Fréchet-derivative of F at (x, y) with respect to x. Moreover, if F is twice Fréchet-differentiable at (x, y), we define In the following, we introduce some concepts of the convergence of set sequences and mapping sequences from [20] which will be used in the next section. Define where N denotes the set of all positive integer numbers.
The continuity properties of a set-valued mapping S can be developed by the convergence of sets.
Definition 2.3. Consider now a family of functions f ν : R n →R, whereR = R ∪ {±∞}. One says that f ν epi-converges to a function f : R n →R as ν → ∞, if the sequence of sets epi f ν converges to epi f in R n × R as ν → ∞.
The characterization of the epi-convergence can be described by the following result.
Next, we recall some basic concepts that are often employed in the literature on the mathematical program with complementarity constraints problem.
Let (x, y) be a feasible point of problem (1.1) and for convenience let us define the index sets The linear independence constraint qualification for SMPCC is as follows.
Definition 2.4. Assume Ψ is continuously differentiable at (x, y). We say the MPCC linear independence constraint qualification (MPCC-LICQ) holds at (x, y) if the vectors of the set are linearly independent, where e i denotes the vector with 1 in the ith component but 0's everywhere else.
For a deterministic MPCC, it is well known that the usual nonlinear programming constraint qualifications such as the Mangasarian-Fromovitz constraint qualification (MFCQ) do not hold (see [28,Proposition 1.1]). Since there are several different approaches to reformulate MPCC, various stationarity concepts arise (see e.g. [21] and [29]). We use the following two stationarity concepts for SMPCC. Definition 2.5. Assume (x, y) is a feasible point of SMPCC (1.1), Ψ(·, ·), E[f (·, ·, ξ)] are continuously differentiable at (x, y). Suppose there exist vectors u ∈ R m and v ∈ R m such that (x, y) satisfies the condition By [21, Lemma 1], the C-stationary condition in Definition 2.5 is the nonsmooth Karash-Kuhn-Tucker (KKT) condition using the Clarke generalized gradient [3] by reformulating SMPCC as a nonsmooth stochastic nonlinear programming problem The S-stationary condition [13] is the KKT condition for the relaxed SMPCC The following upper level strict complementarity condition was used in [21] in the context of sensitivity analysis for MPCC. Definition 2.6. We say that the upper level strict complementarity condition (ULSC) holds at (x, y) ifū i and v i , the multipliers corresponding toȳ i and It is well known that a point (x, y) satisfies the lower level strict complementarity condition (LLSC) if y i + Ψ i (x, y) > 0 holds for all i ∈ {1, . . . , m}. We can see from an example from [21] that the ULSC condition is considerably weaker than the LLSC condition, and in practice, it may make more sense than the latter.
We use the following second order condition based on the MPCC-Lagrangian: of (1.1).

Formulating the smoothing SAA method
The log-exponential function is a smoothing function for max-type functions. Let which was studied by many authors, see Rockafellar and Wets [19], [20], Li [7], [8].
We know from the definition that w(t, x) is a smoothing function with respect to x for t > 0 and hence utilizing this property, over the past decade, some authors have used the log-exponential function to propose smoothing methods for generalized linear complementarity problems, nonlinear complementarity problems, and mathematical programs with complementarity constraints, see [14], [16], [30] and references therein. Notice that G(x, y) = min{Ψ(x, y), y} = 0 can be approximated by the equation and G 0 (x, y) = G(x, y). By taking independently and identically distributed random samples ξ i , i = 1, . . . , N , and introducing the smoothing function G t (·, ·) (3.3), we obtain the approximation of problem (1.1)

Convergence of the smoothing SAA method
Throughout the paper, we assume the sample ξ 1 , . . . , ξ N of the random vector ξ is iid and introduce the following assumptions to make (1.1) more clearly defined and to facilitate the analysis.
Assumption 2. For any (x, y) ∈ R n+m , there exist a neighborhood D of (x, y) and a nonnegative measurable function g(ξ) such that E[g(ξ)] < ∞ and for all ξ ∈ Ξ.
Assumptions 1-2 are popularly used conditions for the analysis of SAA method for stochastic programming. Under these two assumptions, we know from [24, The following lemma results straightforwardly from the Uniform Laws of Large Numbers in [24, Theorem 6.36]. as N tends to infinity, then we obtain

Almost sure convergence of optimal solutions
In this subsection, by using the notion of epi-convergence in [20], we establish the almost convergence of optimal solutions of smoothed SAA problem (3.4) to those of SMPCC (1.1) as the sample size tends to infinity.
Let us introduce some notions: For any positive numbers γ and ε, let N ∈ N be such that √ mt N ln 2 < ε for all N N . Then for all N N we obtain The conclusion follows. Now we give a conclusion about the almost sure convergence of the set Z N as N tends to infinity in the following proposition.
This is almost surely a continuous mapping from the compact convex set B(0, δ) to itself. By Brouwer's fixed point theorem, ϕ has a fixed point w.p. 1. Hence, there exists a vector (x N , y N ) ∈ B((x, y), ε) w.p. 1 such that (x N , y N ) = ϕ(x N , y N ) = z(H N (x N , y N )) w.p. 1. Therefore, we have from (4.2) that w.p. 1, We know from Lemma 4.1 that Then also by Proposition 2.1, we obtain the conclusion. R e m a r k 4.1. From the above proof, we know that if the condition κ 0 being finite is replaced by f being proper and −∞ < inf f < ∞, the conclusion in Theorem 4.1 can also be obtained.

Almost sure convergence of stationary points
In practice, finding a global minimizer might be difficult and in some cases we might just find a stationary point. As a result, we want to know whether or not a cluster point of the sequence of stationary points is almost surely a kind of stationary point of SMPCC (1.1). For this purpose, we need to investigate the almost sure convergence of stationary points of smoothed SAA problem (3.4) with the sample size tends to infinity.
Notice that (3.4) is a standard nonlinear programming with smooth constraints. If (x N , y N ) is a local optimal solution of the smoothed SAA problem (3.4), then under some constraint qualifications, (x N , y N , λ N ) is a stationary point of (3.4) almost surely, namely, there exists a Lagrange multiplier λ N ∈ R m such that the vector (x N , y N , λ N ) satisfies the following Karash-Kuhn-Tucker (KKT) condition for problem (3.4): Inequality   ∇ The following lemma is important for deriving the convergence of our smoothing SAA method for SMPCC.
Notice that by Lemma 4.1, Then we can obtain the conclusions (ii) and (iii) in the same way as in the proof of (i). We have completed the proof.
By Lemma 4.1, we easily get the relationship of MPCC-LICQ between the SAA problem (1.2) and the true problem (1.1) when N is sufficiently large.
are linearly independent almost surely.
We now prove the almost sure convergence of the smoothing SAA method for SMPCC (1.1).
For i ∈ β, due to the boundedness of η N 1i and η N 2i , for simplicity, we assume lim N →∞ η N 1i = η 1i w.p. 1 and lim N →∞ η N 2i = η 2i w.p. 1. Then taking the limit in the equation (4.10) and by Lemma 4.1, we have we obtain the conclusion of (i) from Definition 2.5. Next we prove that under conditions (ii),ū i 0, v i 0 w.p. 1 for i ∈ β. We assume by contradiction thatū j < 0 w.p. 1 for some j ∈ β. Since MPCC-LICQ holds at (x, y) w.p. 1 by Lemma 4.6, we can choose a vector d N ∈ R n+m such that w.p. 1 for N large enough, Then we obtain w.p. 1, We know from (4.8) in Lemma 4.4 that w.p. 1, which, together with (4.11) and by Lemma 4.4, means that for i = j, Similarly, for i = j, we have As a result, combining (4.12), (4.14), and (4.15), we obtain w.p. 1, Since the ULSC condition holds at (x, y) w.p. 1, we have lim Next we show that {d N } can be chosen bounded w.p. 1 for d N satisfying equations (4.11) w.p. 1 for each N . Let H(z) = T (z) − p and T (z) = Az − b, where Then we know from Lemma 4.6 that A has full row rank and hence there existsd such that H(d, 0) = 0. By the implicit function theorem, there exist positive numbers ε, δ and a continuous function z(·) : εB → B(d, δ) such that z(0) =d and H(z(p), p) = 0 for p ∈ εB, which means that T (·) is so called subinvertible [5] at (d, 0). Moreover, We have from Lemma 4.1 that Q N (z) → 0 w.p. 1 as N → ∞, which implies that for a bounded neighborhood U ofd and any positive number σ, when N is large enough, which, together with the subinvertibility of T (·), by [5, Proposition 3.1] means that when N is sufficiently large, Notice that z ∈ J(Q N ) w.p. 1 means A N z = b N w.p. 1. This, together with (4.19), implies that {d N } can be chosen almost surely bounded. By Lemma 4.1, we obtain and for i = 1, . . . , m, This, together with the almost sure boundedness of {d N }, {η N 1i }, and {η N 2i }, i = 1, . . . , m, leads to the almost sure boundedness of As a result, combining (4.16) and (4.17), we can choose a sequence N k ⊂ N such that This contradicts the condition that (x N , y N , λ N ) satisfies the second order necessary conditions almost surely. Hence, we have thatū i 0 w.p. 1 holds for all i ∈ β. Similarly, v i 0 w.p. 1 for all i ∈ β. Therefore, we know from Definition 2.5 that (x, y) is w.p. 1 an S-stationary point of SMPCC (1.1).
Utilizing the MPCC-SSOSC (Definition 2.7), we obtain the following theorem through an application of standard NLP stability theory.
Notice that the equation (4.20) can be seen as the KKT condition of the NLP problem The MPCC-SSOSC ensures the strong second order sufficient condition for the NLP problem (4.22), which, under MPCC-LICQ, implies the stability of (4.22) in the sense of Robinson [17]. Hence, there exist positive numbers ε, δ, c such that for every p ∈ B(0, ε), the mapping Σ(p) = {z ∈ R n+m+|α|+2|β|+|γ| : 0 ∈ Φ(z) + p, z = (x, y, u, v)} has only one solution z(p) ∈ B(z, δ), z = (x, y,ū, v) = z(0) and the mapping z(·) : B(0, ε) → B(z, δ) satisfies We know from the proof of Lemma 4.5 that when δ is sufficiently small then for i ∈ α, As a result, combining (4.25)-(4.29), we obtain that when δ is sufficiently small, which implies that for ε > 0, when N is sufficiently large, Q N 1 δ < ε w.p. 1. In addition, we know from the definition and the Uniform Laws of Large Numbers that when δ is sufficiently small, then sup z∈B(z,δ) which implies that for the above ε > 0, when N is sufficiently large, Note that t N ln 2 → 0 w.p. 1 as N → ∞. Hence, we know from (4.24), (4.30), (4.31), and (4.32) that for the above ε > 0, when δ is sufficiently small and N sufficiently large, then Applying Brouwer's fixed point theorem to the mapping z(Q N (·)) : B(z, δ) → B(z, δ), where z(·) is defined as in (4.23), we conclude that there is at least one fixed point Therefore, when N is sufficiently large, there exists z N ∈ B(z, δ) w.p. 1 such that 0 ∈ Φ(z N ) + Q N (z N ) w.p. 1, namely, then we have from (4.34) that (x N , y N ) is almost surely a stationary point of (3.4) and λ N is the corresponding multiplier. Furthermore, by (4.35), we have (x N , y N ) → (x, y) w.p. 1 as N → ∞. The proof is completed.

Numerical results
In this section we present some preliminary numerical results obtained by the smoothing SAA method based on the log-exponential function (3.1). Our numerical experiments have been carried out in Matlab 7.1 running on a PC with Intel Pentium M of 1.60 GHz CPU and our tests are focused on different values of the smoothing parameter t and sample size N .
In our experiments, we set the initial values of N k and t k as N 1 = 100 and t 1 = 5, respectively. Then we employed the random number generator unifrnd in Matlab 7.1 to generate independently and identically distributed random samples {ξ 1 , ξ 2 , . . . , ξ N k }. We solved problem (3.4) with N = N k and t = t k by the solver fmincon in Matlab 7.1 to obtain the approximate optimal solution (x N k , y N k ). The initial point was (1, . . . , 1) T . In addition, the parameters were updated by N k+1 := min{10N k , 10 5 } and t k+1 := max{0.1t k , 0.005}. Throughout the tests, we recorded number of iterations of fmincon (Iter), the values of the objective function of problem (3.4) at (x N k , y N k ) (Obj) and these quantities are displayed in the tables of test results.

Conclusion and further remarks
In this paper, we propose a smoothing SAA method for a SMPCC by using the log-exponential function. Utilizing the notion of epi-convergence in variational analysis, we establish the almost sure convergence of optimal solutions generated by the smoothed SAA problem. Moreover, under suitable conditions, we show that any cluster point of the KKT point sequence generated from the smoothed SAA problem is almost surely an S-stationary point of SMPCC as the sample size tends to infinity. The preliminary numerical results indicate that the proposed method is able to solve SMPCC successfully.
However, the smoothing SAA method based on the Fischer-Burmeister function [4] in [22] is difficult to be extended for solving this kind of problems.
Since there is no assumption on measurability of the selections considered, the almost sure convergence appears relatively weak. The uniform almost sure convergence [25] is a more convenient convergence in non-measurable case. Whether the almost sure convergence results in this paper can be extended to uniform almost sure convergence results is one of the important topics in our further study.