A SEQUENTIAL CONVEX PROGRAM METHOD TO DC PROGRAM WITH JOINT CHANCE CONSTRAINTS

. In this paper, we consider a DC (diﬀerence of convex) programming problem with joint chance constraints (JCCDCP). We propose a DC function to approximate the constrained function and a corresponding DC program (P ε ) to approximate the JCCDCP. Under some mild assumptions, we show that the solution of Problem (P ε ) converges to the solution of JC-CDCP when ε ↓ 0. A sequential convex program method is constructed to solve the Problem (P ε ). At each iteration a convex program is solved based on the Monte Carlo method, and the generated optimal sequence is proved to converge to the stationary point of Problem (P ε ).

1. Introduction. In this paper, we consider the following optimization problem: where X ∈ R n is a compact convex set, ξ ∈ R k is a parameter vector, f i : R n → R, i = 1, 2 and c j : R n+k → R, j = 1, . . . , m are real-valued convex functions. Let f (x) := f 1 (x) − f 2 (x), the function f (x) is called DC (difference of convex) function, and Problem (1) is a DC programming problem. The DC program (1) has received a lot of interests, since many complicated optimization problems can be 734 XIANTAO XIAO, JIAN GU, LIWEI ZHANG AND SHAOWU ZHANG transformed into DC program. For instance, optimization problem with complementarity constraints, quadratic programming with quadratic constraints [1], low rank optimization problem [4], etc. See [7] for a good survey on applications, mathematical properties and numerical methods of DC programming.
In this paper, we assume that the parameter ξ of Problem (1) is a random vector (which is very common in practical applications), then the problem (1) is a stochastic DC programming problem, which seems impossible to be solved, even a feasible point is very difficult to obtain. In order to overcome the difficulty, we relax the constraints and consider the following problem x ∈ X .
Usually, the function p(x) is not convex and has no closed form. To solve the joint chance constrained problem, many conservative approximations of p(x) have been proposed. Among all convex approximations, the conditional value-at-risk (CVaR) approximation [9] is proved to be the best, which uses to be a conservative approximation of constraint (2), where The constraint CVaR α (c(x, ξ)) ≤ 0 can be rewritten as It is shown that the CVaR approximation is actually using 1 t [t + z] + to approximate 1 (0,+∞) (z). From Figure (1) we can see that, although CVaR approximation is the best convex approximation of p(x), however, it is not a good choice since the approximation becomes terrible as z → +∞. Hong, Yang and Zhang [6] proposed a better but nonconvex approximation, they used a DC function 1 Figure 1. CVaR and DC approximation to approximate 1 (0,+∞) (z). From Figure (1), when t > 0 is very small, the approximation is much better than CVaR approximation. But the drawback is also obvious, it is not convex. By using the DC approximation function to approximate p(x), the JCCDCP can be approximated by a constrained problem with both object function and constraint function are DC functions. In this paper, we show that, under some mild assumptions, the optimal solution of the DC approximation problem is an approximate solution of JCCDCP, when the parameter t is close to 0. To solve the DC programming problem, we propose a sequential convex program method.
The study in this paper gives a deep insight into the work of Hong, Yang and Zhang [6] and extends the results in [6] to a more general form. Compared to the work in the literature, the problem we consider is a stochastic DC programming problem which is more general and more complicated. The assumptions used in [6] are very strict, especially for Assumption 3 "c(x, ξ) is differentiable with respect to x w.p.1", which in fact reduces the joint chance constraint to a single chance constraint. In this paper, we give a more relaxed version than Assumption 3, and remove some other assumptions.
The rest of this paper is organized as follows. In Section 2, we review the DC approximation of the joint chance constraints, and develop a good approximation of JCCDCP, denoted by (P ε ). In Section 3, we propose a sequential convex program method to solve Problem (P ε ). The numerical results are reported in Section 4, followed by the conclusions in Section 5.
2. Approximation of JCCDCP. In this section, we propose an approximation of Problem (P) in the form of (P) which is equivalent to JCCDCP.
Before the discussion, we give the following assumptions.
Then Ω 0 = cl Ω I 0 . Assumption 2.4 is a kind of constraint qualification. It is implied by the Slater condition when Ω 0 is convex.
Under the above assumptions, we have the following results.
Lemma 2.6. If Assumption 2.3 holds, then for any x ∈ X ,p(x, t) is non-decreasing in t when t > 0. Lemma 2.7. Suppose that Assumptions 2.1 to 2.3 are satisfied. Then, for any t ∈ (−δ, δ) and x 0 ∈ O, (i): Proof. Note g 1 (x, t) is convex in (x, t) and x 0 ∈ int (dom g 1 (·, t)), then by [12,Theorem 7.47] and the formula on subgradient of convex function in [10], the result in (i) is easy to be obtained. Let ψ(t) = g 1 (x, t). Then by [12,Theorem 7.44 The following theorem is critical which is the base of εapproximation theory.
Theorem 2.8. Suppose that Assumptions 2.1 to 2.3 are satisfied. Then, (3) Proof. Since g 1 (x, t) is convex and directional differentiable, then Then by Lemma 2.6 and 2.7(ii), The proof is complete.
the approximation problem based on the constraint, denoted by (P ε ), is in the form of x ∈ X .
In the rest of this section, we show that, when ε > 0 is small enough, the DC optimization Problem (P ε ) is a good approximation of Problem (P).
Let Ω(ε) = {x ∈ X : denote the deviation of set A from set B. Let S(ε) and ν(ε) be the set of optimal solutions and the optimal value of Problem (P ε ), S 0 and ν 0 be the set of optimal solutions and the optimal value of Problem (P). Then we have the following results. Theorem 2.9 shows that when ε > 0 is small enough, the optimal solution of Problem (P ε ) is a good approximation to Problem (P). However, since Problem P ε is a DC program, we can only obtain a stationary point. Thus, we should discuss whether the stationary point of Problem (P ε ) converges to the stationary point of (P) when ε → 0.

XIANTAO XIAO, JIAN GU, LIWEI ZHANG AND SHAOWU ZHANG
Let Λ 0 and Λ ε be the sets of stationary points of Problem (P) and (P ε ), respectively, namely, Theorem 2.10. Suppose that Assumptions 2.1 to 2.5 are satisfied. Then According to the outer semicontinuity of N X and∂f , and the definition of Λ ε k , take k → ∞, we have By the definition of Λ 0 we get (x, λ) ∈ Λ 0 . Thus lim sup ε↓0 Λ(ε) ⊂ Λ 0 . The proof is complete.
3. Sequential convex program method for problem (P ε ). Problem (P ε ) is a DC program, because the object function f 1 (x) − f 2 (x) and the left hand side of the constraint g 1 (x, ε) − g 2 (x) ≤ ε(1 − α) are both DC functions. In this section, we propose a sequential convex program method to obtain a stationary point of where w ∈ ∂g 2 (y). Since g 2 (x) is convex, for any y ∈ X and w ∈ ∂g 2 (y), we have Then Ω(ε) y,w ⊂ Ω(ε). (6) For y ∈ X , Ω(ε) y,w is a convex subset of Ω(ε), since Then for y ∈ X, Problem CP(ε, y, w, v) is a conservative convex approximation of Problem (P ε ). By solving a sequence of Problem CP(ε, y, w, v), we can obtain a stationary point of Problem (P ε ).
Step 1: where N R− (g(x k , ε)) is the norm cone of R − at g(x k , ε), and Step 2: Choose v k ∈ ∂f 2 (x k ),w k ∈ ∂g 2 (x k ), solve Problem CP(ε, x k , w k , v k ) to obtain its optimal solution x k+1 .

Remark 1.
For Algorithm 1, the following works in the literature should be mentioned.
• For the case f 2 = 0, i.e., the object function is convex, Hong et al. [6] provided a comprehensive analysis on the convergence of the algorithm. • In [3], a low rank constrained optimization problem is discussed, and the technique by linearizing the second convex term of DC function also plays an important role in the algorithm design.
To get the convergence result of the above algorithm, we need the following assumption.

XIANTAO XIAO, JIAN GU, LIWEI ZHANG AND SHAOWU ZHANG
Define the indicator function δ X (x) as Lemma 3.3. Suppose that Assumptions 2.1 to 2.4 are satisfied, and Assumption 3.2 holds. If g 2 is strictly convex on X , then for ε > 0 small enough, y ∈ Ω(ε) and w ∈ ∂g 2 (y), there exists γ > 0 such that where Since X is compact, we assume that y i → y, w i → w and x i → x. Next we prove y = x. Because x i ∈ bdΩ(ε i ) y i ,w i , we can assume that Assume that ξ i → ξ, η i → η, then ξ ∈ lim sup t↓0,x →x ∂ x g 1 (x , t), η ∈ N X ( x), w ∈ ∂g 2 ( x).
Then we have the following theorem on the convergence of the sequential convex program method. Proof. Suppose that (x, w, v) is a cluster point of the sequence {x k , w k , v k }. By Lemma 3.1, we obtain {x k } ⊂ Ω(ε), then the cluster point of {x k } satisfiesx ∈ Ω(ε). Since ∂f 2 and ∂g 2 is outer semicontinuous atx, there exists δ 1 > 0, such that for j ∈ N large enough, x j ∈ {x k } ∩ B δ1 (x), ∂f 2 (x j ) and ∂h(x j ) are bounded, v ∈ ∂f 2 (x) and w ∈ ∂g 2 (x).
By Lemma 3.1, {f (x k )} is monotone non-increasing, which converges to inf k f (x k ), then f (x) = inf k f (x k ). Since the Assumption 3.2 is satisfied, the condition (8) holds for Ω(ε)ȳ ,w . We first provex is an optimal point of Problem CP(ε,x, w, v). Let {k j } is an infinite sequence of N such that (x kj , w kj , v kj ) → (x, w, v). Definef By [11,Proposition 7.4], δ Ω(ε) x k j ,w k j (·) epi-converges to δ Ω(ε)x,w (·), then we get f kj (·) epi-convergesf (·). By [11,Theorem 7.33], we obtain Note that argminf kj is the optimal set of Problem CP(ε, x kj , w kj , v kj ), argmin f is the optimal set of Problem CP(ε,x, w, v). By (17), the sequence {x kj +1 } has a cluster pointx ∈ Ω(ε), which is in the set argminf . Since f (x) ≤ f (x) and . Thenx is the optimal solution of Problem CP(ε,x, w, v). Because Problem CP(ε,x, w, v) is a convex problem and Slater condition holds, we have which impliesx is a stationary point of Problem (P ε ).

Numerical experiments.
In this section, we will first consider some numerical issues for implementing our algorithm and then show the results for some numerical examples.

Monte Carlo method.
To implement Algorithm 1, at each iteration, we need to solve a convex program The problem is difficult to solve because g 1 (x, ε) and g 2 (x) are expressed by mathematical expectations, and the closed forms of them are unknown. To overcome this difficulty, we apply a Monte Carlo method. Recall the definitions of g 1 , g 2 and Lemma 2.7, we have Let ξ 1 , . . . , ξ n denote an independent and identically distributed (i.i.d.) sample of ξ. Then, a natural estimator of g(x) is We consider the problem and use its optimal solution to approximate the optimal solution of Problem (19). Let S and ν * denote the set of optimal solutions and the optimal objective value of Problem (19), andŜ n andν * n denote the set of optimal solutions and the optimal objective value of Problem (20). Then we have the following lemma to show that S n andν * n converge to S and ν * . Lemma 4.1. [6, Theorem 6] Suppose that Assumptions 2.1 to 2.5 are satisfied. When ε > 0 is small enough, D(Ŝ n , S) → 0 w.p.1 andν * n → ν * w.p.1 as n → ∞. 4.2. Numerical environment. The numerical experiment is carried out in MAT-LAB 7.10 running on a PC Intel Core i5 CPU (3.20 GHz) and 4 GB RAM. In the following test problems, we set α = 0.9, ε = 0.05 2 . Let {x k } be the generated sequence. We stop the algorithm if norm(x k+1 − x k ) 2 ≤ 10 −6 .

4.3.
A L 1 optimization problem. We consider the following problem: where ξ ij , i = 1, . . . , m and j = 1, . . . , n are independent and identically distributed standard norm random variables, b is a constant. The Problem (21) can be written equivalently as in the sense that Problem (21) and Problem (22) have the same optimal value, and the optimal set for Problem (22) is a subset of that for Problem (21). Note that Problem (22) is a simple type of JCCDCP (the object function is linear).
Lemma 4.2. The optimal solution x * of Problem (22) is is the inverse chi-squared distribution function with n degrees of freedom.
Proof. By the symmetric form of Problem (22), the optimal solution x * satisfies x * 1 = x * 2 = · · · = x * n . Then the joint chance constraint of Problem (22) in x * takes the form   Pr which is equivalent to Then, it is clear that We implement Algorithm 3 for the following two problems. At each iteration in the algorithm, we use CVX [5] to solve a convex program generated by the Monte Carlo method. Problem 4.3. In Problem (22), let n = 3, m = 2, b = 9, N = 100. By Lemma 4.2, the optimal solution is (1.08, 1.08, 1.08) T , the optimal value is -3.23.
The numerical results are shown in Figure 2. Take Problem 4.4 for example, the algorithm typically requires less than 20 iterations to converge to the optimal value. When sample size N = 100, at each iteration, a convex program with 4223 variables and 3112 constraints is solved in 26 seconds on average. When N grows, the size of the convex program increases much faster, and we found CVX performed not well under this situation.

4.4.
A nonconvex quadratic programming problem. Let G ∈ R n×n be a symmetric but not positive semidefinite matrix. We consider the following nonconvex stochastic quadratic programming problem, where with ξ i = (ξ i1 , . . . , ξ in ) be a m × n matrix of random variables, c, l, u ∈ R n and b ∈ R m . Assume that λ is the largest eigenvalue of G. Let ρ > max(0, λ) be a positive number and I ∈ R n×n be the identity matrix. Define Note that both f 1 and f 2 are convex. Then the Problem (24) can be reformulate as which is a JCCDCP.
The numerical results for Problem 4.5 are reported in Table 1, where "It.", "Res0." and "Res*." stand for the number of iterations, the residuals φ(·, ·, ·) at the first iteration and at the final iteration, respectively. For Problem 4.5, we use fmincon, which is a function from MATLAB Optimization Toolbox, to solve the convex program at each iteration in the algorithm. Problem 4.5 is a nonconvex stochastic program and the optimal solution is no longer known. However, from the numerical results in Table 1, by our method we can always obtain a stationary point in a suitable time even for a medium scale problem. Also from Table 1, as the sample size N grows, the size of the convex program in each iteration increases fast, consequently the CPU time becomes much longer.

5.
Conclusions. In this paper, we introduced a DC technique to approximate the DC programming problem with joint chance constraints. Under some assumptions, we had shown the stationary point of the approximation problem converges to a stationary point of the original problem when the parameter ε ↓ 0. We proposed a sequential convex program method for solving the approximation problem and analyzed its convergence. Numerical experiments conducted on a variety of JCCDCP problems demonstrate that our method is efficient.
At each iteration, it needs to solve a large scale convex program because we use the Monte Carlo method. From our experience, neither CVX nor fmincon is a satisfactory choice to solve this type of problem rapidly and efficiently. Thus, how to solve the large scale convex program rapidly is a key factor that influences the efficiency of our method. On the other hand, many robust optimization algorithms have been proposed to avoid using the Monte Carlo method. Whether these robust optimization techniques can be applied in our method is also attractive.