Solving Chance-Constrained Problems via a Smooth Sample-Based Nonlinear Approximation

TL;DR

This paper presents a novel approach for solving chance-constrained optimization problems by reformulating probabilistic constraints as quantile functions and approximating them with a differentiable sample average, enabling the use of standard solvers.

Contribution

The authors introduce a differentiable sample-based approximation of the quantile function for chance constraints, with theoretical guarantees and a new trust-region method for joint constraints.

Findings

Method scales well with sample size.

Smoothing reduces bias in chance constraint approximation.

Empirical results demonstrate effectiveness on various problems.

Abstract

We introduce a new method for solving nonlinear continuous optimization problems with chance constraints. Our method is based on a reformulation of the probabilistic constraint as a quantile function. The quantile function is approximated via a differentiable sample average approximation. We provide theoretical statistical guarantees of the approximation, and illustrate empirically that the reformulation can be directly used by standard nonlinear optimization solvers in the case of single chance constraints. Furthermore, we propose an S$\ell_1$QP-type trust-region method to solve instances with joint chance constraints. We demonstrate the performance of the method on several problems, and show that it scales well with the sample size and that the smoothing can be used to counteract the bias in the chance constraint approximation induced by the sample approximation.

Tables

Table 1. Table 1 : Iterations performed by Knitro when maximizing x 𝑥 x using two different formulations to describe the feasible region and different bounds on x 𝑥 x . ∗ : Knitro gives the following error message: “Convergence to an infeasible point. Problem appears to be locally infeasible.”

Bounds	$[-1,1]$	$[-10,10]$	$[-100,100]$	$(-\infty ,\infty )$
$\mathbb{P}\left(C(x,\xi )\le 0\right)\ge 1-\alpha$	6	36	14	$5^{\ast }$
$Q\left(C(x,\xi )\right)\le 0$	6	6	6	6

Table 2. Table 2 : Number of iterations and CPU time needed by Knitro to return a solution.

N	Avg. time (sec)	Max. time (sec)	Avg. iter (per $ϵ$)	Avg. $ϵ$
200	1.1301	1.8285	13.236	1628.3
500	1.4797	2.6307	14.250	1135.6
2000	10.293	13.643	13.479	700.68
5000	49.943	113.91	14.073	554.47

Table 3. Table 3 : Statistics for the multi-knapsack problem.

Sample size	100		500		1,000
Method	NLP	MIP	NLP	MIP	NLP	MIP
Min. prob	0.9499	0.9445	0.9499	0.9503	0.9499	0.9502
Avg. prob	0.9504	0.9591	0.9503	0.9520	0.9511	0.9514
Max. prob	0.9521	0.9723	0.9526	0.9549	0.9591	0.9535
Min. obj	5816.6	5771.3	5834.4	5830.3	5827.5	5831.3
Avg. obj	5828.9	5806.0	5839.0	5833.2	5840.2	5836.6
Max. obj	5842.6	5836.1	5842.6	5839.2	5842.7	5841.0
Avg. time (s)	8.7	90	29.3	168	37.4	372
Max. time (s)	24.6	120	49.8	180	80.5	420
Avg. $ϵ$	17	-	9.8	-	7.5	-

Equations

$$\underset{x\in X}{\min }$$

$$s.t.$$

$$Q^{1-\alpha }(Y)=\inf \{y\in \mathbb{R}\mid P(Y\le y)\ge 1-\alpha \}.$$

$$\underset{x\in X}{\min }$$

$$s.t.$$

$${\tilde{F}}^N(t;x)=\frac{1}{N}\sum _{i=1}^N\text{mathbbm}1(C(x,{\xi }_i)\le t).$$

$${\tilde{Q}}^{1-\alpha }(C^N(x))=\inf \left\{y\mid \frac{1}{N}\sum _{i=1}^N\text{mathbbm}1\left(C(x,{\xi }_i)\le y\right)\ge 1-\alpha \right\}=C_{[M]}(x),$$

$$\left\{x\in {\mathbb{R}}^2\mid Q^{1-\alpha }(C(x,\xi ))\le 0\right\},$$

$$F_{ϵ}^N(t;x)=\frac{1}{N}\sum _{i=1}^N{\Gamma }_{ϵ}(C(x,{\xi }_i)-t),$$

$$\displaystyle\Gamma_{\epsilon}(y)=\left\{\begin{array}[]{lr}1,&y\leq-\epsilon\\ \gamma_{\epsilon}(y),&-\epsilon<y<\epsilon\\ 0,&y\geq\epsilon\end{array}\right.$$

$$\sum _{i=1}^N{\Gamma }_{ϵ}([C^N(x){]}_i-Q_{ϵ}^{1-\alpha })=(1-\alpha )N,$$

$$\underset{x\in X}{\min }$$

$$s.t.$$

$${\hat{F}}_{ϵ}^N(t;x)$$

$$\mathbb{E}\left[{\left({\hat{F}}_{ϵ}^N(t;x)-F(t;x)\right)}^2\right]\approx \frac{F(t;x)\left[1-F(t;x)\right]}{N}-u(t)+{ϵ}^{4}v(t),$$

$$u(t)=f(t;x)\frac{ϵ}{N}\left(1-{\int }_{-1}^1{\Gamma }_1^2(y)dy\right),v(t)={\left(\frac{1}{2}{f}^{\prime }(t;x){\int }_{-1}^1{y}^2{\Gamma }_1^{\prime }(y)dy\right)}^2$$

$${\gamma }_{ϵ}(y)=\frac{15}{16}\left(-\frac{1}{5}{\left(\frac{y}{ϵ}\right)}^5+\frac{2}{3}{\left(\frac{y}{ϵ}\right)}^3-\left(\frac{y}{ϵ}\right)+\frac{8}{15}\right),$$

$$X_{\alpha }=\left\{x\in X\mid F(0;x)\ge 1-\alpha \right\},$$

$$F_{ϵ}(t;x)={\int }_{{\mathbb{R}}^s}{\Gamma }_{ϵ}\left(C(x,\xi )-t\right)d{P}_{\xi },$$

$$X_{ϵ,\delta }^{N,t}=\left\{x\in X\mid F_{ϵ}^N(-t;x)\ge 1-\delta \right\},\delta \in [0,1],t\in \mathbb{R}.$$

$$\underset{x\in U}{\sup }\left|F_{ϵ}^N(0;x)-F(0;x)\right|\to 0,\text{w.p.1},\text{as }N\to \infty \text{ and }ϵ\to 0.$$

$$\underset{x\in U}{\sup }\left|F_{ϵ}^N(0;x)-F(0;x)\right|<\delta .$$

$$\underset{k\to \infty }{\mathrm{liminf}}{v}_{{ϵ}_k}^{N_k}\ge v^{\ast },\text{w.p.1}.$$

$$\underset{k\to \infty }{\mathrm{limsup}}{v}_{{ϵ}_k}^{N_k}\le f(\bar{x})=v^{\ast },\text{w.p.1}.$$

$$P\left\{\sum _{1=1}^N\left(Y_i-E[Y_i]\right)\ge tN\right\}\le \exp \left\{-\frac{2N^2{t}^2}{{\sum }_{i=1}^N(b_i-a_i{)}^2}\right\}.$$

$$P\left(F_{ϵ}^N(-t;x)\ge 1-\delta \right)\le \exp \left\{-2N{M}_x^2\right\}.$$

$$P\left(F_{ϵ}^N(-t;x)\ge 1-\delta \right)$$

$$\le P\left(F_{ϵ}^N(-t;x)-F_{ϵ}(-t;x)>F(0;x)-F_{ϵ}(-t;x)+(\alpha -\delta )\right)$$

$$=P\left(F_{ϵ}^N(-t;x)-F_{ϵ}(-t;x)>M_x\right)$$

$$=P\left(\sum _{i=1}^N\left(Y_i-E[Y_i]\right)\ge M_{x}N\right)$$

$$\le \exp \left\{-\frac{2N^2{M}_x^2}{{\sum }_{i=1}^N{\left(1-0\right)}^2}\right\}=\exp \left\{-2N{M}_x^2\right\},$$

$$P\left(F_{ϵ}^N(0;x)\ge 1-\alpha \right)\le \exp \left\{-2N{\beta }_x^2\right\},$$

$$F(0;x)-F_{ϵ}(0;x)$$

$$={\int }_{-\infty }^{\infty }\text{mathbbm}1\left(y\right)h_x(y)dy-{\int }_{-\infty }^{\infty }{\Gamma }_{ϵ}(y)h_x(y)dy$$

Full text

{AMS}

90C15, 90C30, 90C55

1 Introduction

Consider the following minimization problem subject to a chance constraint

[TABLE]

where $X\subseteq\mathbb{R}^{n}$ represents a closed deterministic feasible region, $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ represents the objective function to be minimized, $c:\mathbb{R}^{n}\times\mathbb{R}^{s}\rightarrow\mathbb{R}^{m}$ is a vector-valued function, $\xi$ is a random vector with support $\Xi\subseteq\mathbb{R}^{s}$, and $\alpha\in(0,1)$ is a risk parameter given a priori. This kind of problem is known in the literature as a joint chance-constrained problem (JCCP), in which all $m$ constraints must be satisfied simultaneously with probability at least $1-\alpha$. A special case is when $m=1$ and only one constraint needs to be satisfied with probability $1-\alpha$; this is called a single chance-constrained problem (SCCP).

Chance-constrained optimization problems were introduced in [9] and have been extensively studied ever since; see e.g., [37]. These problems can be found in a plethora of applications such as water management [27], optimization of chemical processes [21, 22], and optimal power flow [4, 32], to mention a few.

This formulation gives rise to several difficulties. First, the structural properties of $c(x,\xi)$ may not be passed to the constraint $\operatorname{\mathbb{P}}(c_{j}(x,\xi)\leq 0,j\in\{1,\ldots,m\})\geq 1-\alpha$. For example, even if $c_{1}(x,\xi),\dots,c_{m}(x,\xi)$ are all linear in $x$, the probabilistic constraint may not define a convex feasible region. Second, the functional form of the distribution of $\xi$ may not be known. And finally, even if the distribution of $\xi$ is known, there might not be a tractable analytical function to describe the constraint or a numerically efficient way to compute it.

We focus on problems in which $f$ and $c$ are continuous and differentiable, with respect to $x$, and where $\xi$ is such that the random variable defined by $c(x,\xi)$ has a continuous cumulative density function (cdf) for all $x\in X$. Other than that, we make no assumptions on the type of constraint $c(x,\xi)$; we do not limit ourselves to linear or convex functions.

1.1 Literature Review

In this section we review several methods for solving chance-constrained problems. The first couple of methods make assumptions on the distribution of $\xi$, or at least require that some knowledge of the distribution $\xi$ is available, whereas the latter methods make no such assumptions. Our idea falls into the latter category, and does not impose restrictions on the form of the constraints or on the distribution, except that it is a continuous distribution.

To compute Eq. 1b, some methods take advantage of the special structure of constraint $c$ and the distribution of $\xi$. For example, in [19, Lemma 2.2], if $\xi$ has an elliptical symmetric distribution (e.g., multivariate normal random variable) and if $c(x,\xi)=\xi^{T}x-b$, then the problem can be expressed as a second-order cone program, and can be solved in a tractable and efficient manner using standard solvers. More sophisticated procedures make use of efficient numerical integration techniques to compute Gausian-like probabilities, as well as their corresponding gradients (e.g., [13, 16, 17, 43]). For linear chance constraints, approaches based on cutting planes [45], bundle [46], and SQP [5] methods have been developed. Nonlinear constraints can be addressed using a spherical-radial decomposition method [18, 44].

There are several approaches that avoid making assumptions on $\xi$ by sampling to approximate the probabilistic constraint. We divide such techniques into those that use an integer reformulation, and those that are structured as nonlinear programming problems (NLP).

In the first category, the probabilistic constraint is replaced by the empirical distribution function and formulated as a mixed-integer nonlinear program (MINLP), e.g. [29]. This method takes realizations of the random variable $\xi$ and requires that at least $(1-\alpha)\%$ of the realized constraints are satisfied. This idea has several advantages; the solution converges to the true solution of Eq. 1 as the sample size increases, and for convex instances, a global solution of the approximation can be computed by branch and bound. However, the complexity of the problem increases with the sample size because the number of binary variables increases. This issue has been circumvented to some extent for certain types of constraints $c(x,\xi)$, e.g. linear constraints [30], but this procedure cannot be applied to general constraints.

To avoid the need for solving MINLPs, a conservative and more tractable approximation can be obtained by using the scenario-based approach, e.g. [6, 7, 34]. These methods seek to satisfy all of the constraints in a sample of pre-specified size to ensure that the solution is feasible with high confidence. The reformulation proposed in [6] allows one to obtain solutions very easily since the structure of the original constraints $c(x,\xi)$ is preserved. For example, if $f$ is convex, $X$ defines a convex set and $c(x,\xi_{i})$ define convex functions on $x$ for all $i=\{1,\dots,N\}$ (where $N$ is the size of the sample), then the problem that arises is convex. However, the limitation of this method is that one may obtain highly conservative solutions and that it does not offer direct control on the probability level of the chance constraint.

Another tractable conservative approximation can be obtained by using conditional value-at-risk, which is based on using a conservative convex approximation of the indicator function [33]. This approximation typically finds feasible but sub-optimal solutions. To avoid overly conservative solutions, some authors propose a difference-of-convex functions approximation to the indicator function [25], while some others propose outer or inner approximations to the indicator function [8, 15]. The main drawback of the last two methods mentioned is that the constraint function can become very “flat” in the regions of interest. That is, when the probability of the constraint satisfaction evaluated at a given $x$ is close to one or close to zero, the information obtained from the gradient may not provide a sufficiently accurate approximation of the constraints or may even be erroneous. We discuss this issue further when we motivate our proposed reformulation of Eq. 1 in Section 2.

A potential drawback of sample-based methods is that in some cases the number of samples needed to obtain an accurate solution might be fairly large [20, Example 1]. Tailored approaches that take the problem structure directly into account can be more efficient when applicable.

1.2 Our Approach

In this paper we propose an NLP reformulation of the problem based on sample average approximation (SAA), in which we rewrite the probability constraint as a quantile constraint. We first argue that the latter formulation is more suitable for gradient-based optimization methods than the former. Then, we present a way to approximate the quantile constraint using samples in a way that provides a smooth constraint. We find in our numerical experiments that the smoothing helps control the variance of the solutions obtained from different samples. The approximation we present uses a parameter $\epsilon$ that depends on the scaling of the problem, which might be difficult to select; thus, we also discuss an algorithm that finds an appropriate value of $\epsilon$.

The proposed approach has the following advantages. First, we observe empirically that the solutions we obtain are very often “more robust” compared to solutions obtained from solving the traditional (non-smoothed) SAA problem via mixed-integer programming (MIP) formulations that use the empirical cdf. More specifically, when we compare the solutions obtained from our method to those returned by the MIP, we consistently find that our method attains a better objective function value and has less variability across different samples, while still maintaining out-of-sample feasibility. We attribute this to the fact that our approach uses information from the entire sample, whereas the empirical cdf only takes into account whether enough constraints are satisfied. These numerical experiments suggest that the solutions we attain are able to “generalize”, i.e., the solutions obtain better performance on out-of-sample scenarios. Second, our method is scalable; we find that the time to solve an instance with our approach grows modestly with the sample size. Third, our formulation does not assume any type of structure from the constraint functions save for continuity and differentiability with respect to $x$; therefore, we are able to find local solutions to non-linear and non-convex problems. Finally, for SCCPs, a standard NLP solver can be used to find the solution of the problem for the appropriate choice of the parameters involved in the formulation, and for JCCPs we are able to derive an S$\ell_{1}$QP algorithm that has proven empirically to give fast convergence to the solution after only a few iterations.

The remainder of this paper is organized as follows. In Section 2 we present the proposed reformulation of the problem. In Section 3 we study the convergence and feasibility of the solutions to the approximate problem to those of the true problem Eq. 1. In Section 4 we discuss single chance-constrained problems, and in Section 5 we extend this approach to solve joint chance-constrained problems. In Section 6 we present numerical results. Finally, we make some concluding remarks in Section 7.

2 Reformulation and approximation of the CCP

The quantile of a random variable $Y$ at the level $1-\alpha\in(0,1)$ is defined as

[TABLE]

Therefore, $\operatorname{\mathbb{P}}(C(x,\xi)\leq 0)\geq 1-\alpha$ is equivalent to $Q^{1-\alpha}(C(x,\xi))\leq 0$, where $C:\mathbb{R}^{n}\times\mathbb{R}^{s}\rightarrow\mathbb{R}$ is a real valued function, $x\in\mathbb{R}^{n}$, and $\xi$ is a random vector taking values in $\mathbb{R}^{s}$. Here, we make a distinction between $c(x,\xi)$ and $C(x,\xi)$; $C(x,\xi)$ denotes a scalar-valued function and $c(x,\xi)$ is the vector-valued function defined in problem Eq. 1. The reason for this distinction is that $Q^{1-\alpha}(C(x,\xi))$ is only well defined if $C(x,\xi)$ is a scalar random variable. However, when $c$ is a vector-valued constraint, we set $C(x,\xi)=\max_{j=1,\ldots,m}c_{j}(x,\xi)$ and obtain the following reformulation of problem Eq. 1,

[TABLE]

The advantage of rewriting the constraint as the quantile is that this reduces the “flatness” encountered when the probability function $\operatorname{\mathbb{P}}(C(x,\xi)\leq 0)$ is used as the constraint, see Fig. 1 as an example. Now the feasible region is measured in the image of $C(x,\xi)$ and not in the bounded space $[0,1]$, which is the image of the probability function.

Fig. 1 also shows how the feasible regions coincide for both formulations. Yet, although the two formulations are equivalent, solvers that depend on gradient information are likely to perform better when the quantile function is used, as we will see in the numerical experiments. This can be attributed to the fact that the linearization of the quantile function provides an approximation of the constraint that is reasonably accurate for a wider range than the linearization of the probabilistic constraint. The consequence of “flatness” is that a non-specialized gradient-based solver that starts at a point $x$ where $\operatorname{\mathbb{P}}(C(x,\xi)\leq 0)\approx 1$ takes a big step to decrease the objective, disregarding the constraint function. Thus, the solver is forced to reject many trial steps, which can lead to an increased number of iterations. This observation has been made before in [46, p.734], where the authors state that the chance constraint “can define a gully-shaped feasible region that makes [the chance-constrained problem] hard to tackle with standard nonlinear programming solvers [31].” Moreover, [31] shows that this characteristic becomes more noticeable, at least for normal random variables, as the dimension of the random variable increases.

To illustrate this claim, we compared the numerical performance of the nonlinear programming solver Knitro for each of the two formulations for the constraints $C(x,\xi)=x^{2}-2+\xi$. We used Knitro to maximize $x$ using: (1) the probabilistic constraint $\operatorname{\mathbb{P}}(C(x,\xi)\leq 0)=\Phi(2-x^{2})\geq 1-\alpha$, and (2) the quantile constraint $Q^{1-\alpha}(C(x,\xi))=x^{2}-2+\Phi^{-1}(1-\alpha)\leq 0$; see Fig. 1. We chose $x_{0}=3$ as starting point. Table 1 illustrates the performance in terms of iteration count using the two types of constraints. Each column shows the number of iterations using different bounds on $x$. The results in the table suggest that the quantile function is more robust since the number of iterations is independent of the bounds and the solver always returned feasible and optimal solutions.

Motivated by the aforementioned arguments, we propose to work with the quantile formulation. However, as with the probabilistic formulation, an explicit expression to compute the quantile is usually not available or can be intractable to compute. To address this issue, we introduce a sample average approximation of the quantile in the following section.

Alternative approaches to overcome the convergence failures arising from the flatness of the probability constraints include a logarithmic transformation [31] or the use of interpolation steps [42] within a tailored optimization algorithm.

2.1 Smooth approximation of the quantile function

A quantile estimator can be obtained by inverting an estimator of the cdf when the distribution is continuous. Let $\{\xi_{1},\ldots,\xi_{N}\}$ be an independent and identically distributed (i.i.d.) sample and consider the empirical cdf of $\{C(x,\xi_{1}),C(x,\xi_{2}),\ldots,C(x,\xi_{N})\}$ at a given point $x$,

[TABLE]

An estimate of the $(1-\alpha)$-quantile (and hence $Q^{1-\alpha}(C(x,\xi))$) is obtained from a value $t$ such that $\widetilde{F}^{N}(t;x)\approx 1-\alpha$. The $(1-\alpha)$-empirical quantile at $x$ is defined as

[TABLE]

for $M=\lceil(1-\alpha)N\rceil$. Here, $\mathbbm{1}$ is the indicator function and $C_{[j]}$ denotes the $j$th smallest observation of the values $\{C(x,\xi_{1}),\ldots,C(x,\xi_{N})\}$ for a fixed $x$.

There are two main drawbacks to using the empirical quantile as an approximation. First, the empirical approximation is not differentiable. The “active scenario” (the scenario $\xi_{i}$ such that $C_{[M]}(x)=C(x,\xi_{i})$) changes for different values of $x$. Hence, even if the constraint $C(x,\xi)$ is smooth for a fixed value of $\xi$, this does not hold for $\widetilde{Q}^{1-\alpha}(C^{N}(x))$. Second, the aforementioned feature also introduces artificial local minima to the problem, as can be seen in Fig. 2. For a particular sample $\{\xi_{1},\ldots,\xi_{N}\}$, this figure depicts the feasible region

[TABLE]

and its approximation based on the empirical quantile function; where $C(x,\xi)=\xi^{T}x-1$, $\xi\sim N(\mu,\Sigma)$, $\mu=\mathbf{0}\in\mathbb{R}^{2}$, and $\Sigma$ is the $2\times 2$ identity matrix.

Notice that even if the true feasible region is convex and has a smooth boundary (dotted curve), its approximation (solid curve) is non-smooth and introduces a lot of “inward kinks”, which break the convexity of the feasible region. These “inward kinks” can potentially cause a local optimization algorithm to converge to inferior local minima more easily. Furthermore, as the number of scenarios increases, the number of non-smooth points increases, leading to more sub-optimal local minima.

In order to avoid the aforementioned problems, a smooth approximation of the quantile is obtained by estimating the cdf with a smooth function $\hat{F}$, and solving the equation $\hat{F}(t;x)=1-\alpha$, where $t$ is the approximation of $Q^{1-\alpha}(C(x,\xi))$. Following an approach similar to [15] or [40, Section 4.4.2], we approximate $F$ at a given point $x$ as

[TABLE]

where $\epsilon>0$ is a parameter of the approximation,

[TABLE]

and $\gamma_{\epsilon}:[-\epsilon,\epsilon]\rightarrow[0,1]$ is a symmetric and strictly decreasing function such that it makes $\Gamma_{\epsilon}$ differentiable. With this choice of $\gamma_{\epsilon}$, $F^{N}_{\epsilon}(t;x)$ is a differentiable approximation of the empirical cdf, $\widetilde{F}^{N}(t;x)$.

Assuming that $C(x,\xi)$ is a continuous random variable for a fixed value of $x$, we can approximate the $(1-\alpha)$-quantile as the inverse of $F^{N}_{\epsilon}$ at the $1-\alpha$ level. Thus, we define our $\epsilon$-approximation to the $(1-\alpha)$-quantile at $x$ as the value $Q_{\epsilon}^{1-\alpha}$ such that the following equation holds

[TABLE]

where $C^{N}(x)\in\mathbb{R}^{N}$ is a vector with elements $[C^{N}(x)]_{i}=C(x,\xi_{i})$. However, the value $Q_{\epsilon}^{1-\alpha}$ is not unique if the equation $F^{N}_{\epsilon}(t;x)=1-\alpha$ has more than one solution. In order to guarantee the existence of a unique solution for all possible choices of $\epsilon>0$, we need to ensure that $[F^{N}_{\epsilon}]^{\prime}(t;x)\neq 0$ when evaluated at the quantile of interest. The following proposition gives conditions under which $Q_{\epsilon}^{1-\alpha}$ is uniquely defined.

Lemma 2.1**.**

*Let $\gamma_{\epsilon}:[-\epsilon,\epsilon]\rightarrow[0,1]$ be a strictly decreasing function such that $\Gamma_{\epsilon}$ is differentiable, and let $(1-\alpha)N\notin\mathbb{Z}$. Then, equation Eq. 7 has a unique solution. *

Proof 2.2**.**

Let $\psi(Q)=\sum_{i=1}^{N}\Gamma_{\epsilon}(C(x,\xi_{i})-Q)-(1-\alpha)N$. We have that $\psi^{\prime}(Q)=-\sum_{i=1}^{N}\Gamma^{\prime}_{\epsilon}(C(x,\xi_{i})-Q)$. Given that $\gamma_{\epsilon}$ is strictly decreasing and $\Gamma^{\prime}_{\epsilon}(y)=0$ for all $y\geq\epsilon$ and $y\leq-\epsilon$, we conclude that $\psi^{\prime}(Q)\geq 0$ for all $Q$, which means that $\psi$ is an increasing function. If we prove that $\psi^{\prime}(Q)\neq 0$ for all $Q$ such that $\psi(Q)=0$, then we show that there is a unique solution of equation Eq. 7.

*Assume $\psi^{\prime}(Q)=0$, then $Q\leq C(x,\xi_{i})-\epsilon$ or $Q\geq C(x,\xi_{i})+\epsilon$ for all scenarios $i\in\{1,\ldots,N\}$. This means that $\sum_{i=1}^{N}\Gamma_{\epsilon}(C(x,\xi_{i})-Q)\in\mathbb{N}$. Because $(1-\alpha)N\notin\mathbb{N}$, it means that $\psi(Q)\neq 0$, and therefore any $Q$ such that $\psi^{\prime}(Q)=0$ cannot be a solution. *

Remark 2.3**.**

*The assumption $(1-\alpha)N\notin\mathbb{Z}$ of Lemma 2.1 is mild because if it is not satisfied one can either increase $N$, perturb $\alpha$, or add a constant $b\in(0,1)$ to the left hand side of Eq. 7. In fact, if a constant is added to the left-hand side, the impact will disappear as $N\rightarrow\infty$. *

Lemma 2.1 shows that Eq. 7 indeed implicitly defines a function $Q^{1-\alpha}_{\epsilon}$ that maps a vector $C^{N}(x)\in\mathbb{R}^{N}$ to the root of Eq. 7. In order to find solutions for Eq. 1, we focus on solving the following problem

[TABLE]

The above formulation is a sample average approximation of the original problem Eq. 1. Convergence and feasibility of this approximation are analyzed in Section 3.

One could have approximated constraint Eq. 1b with $F^{N}_{\epsilon}(0;x)\geq 1-\alpha$, where $F^{N}_{\epsilon}$ is given by Eq. 2. However, as mentioned at the beginning of this section and as shown in the numerical results (section 6.2.2), working with the probabilistic constraint, or its estimate, instead of the quantile leads to more iterations or failure of the optimization algorithm. Therefore, in contrast with previous works that express the chance constraint using variations of Eq. 2, we only use $F_{\epsilon}^{N}$ to obtain an approximation of the quantile and work directly with the quantile function.

While other options are possible, our choice of the smoothed indicator function, $\Gamma_{\epsilon}$, is motivated by a kernel estimation of the cdf. In [1, p. 326], the authors present the following estimator

[TABLE]

where $\vartheta(\cdot)$ is a real-valued non-negative function, known as kernel function, that integrates to one over all of its domain (features that make it a probability density function) and is symmetric around 0. The function $\Theta(\cdot)$ is the cdf of $\vartheta(\cdot)$.

If $\vartheta(\cdot)$ has a finite support $[-1,1]$, then $\Theta\left(\frac{t-c(x,\xi)}{\epsilon}\right)$ has the structure of the $\Gamma_{\epsilon}$ function in Eq. 6. Furthermore, according to [2, p. 326], if $\epsilon\rightarrow 0$ as $N\rightarrow\infty$, it can be proven that asymptotically

[TABLE]

where

[TABLE]

The first two terms in Eq. 9 correspond to the asymptotic variance and the third term is the squared asymptotic bias [1, p. 257]. This result indicates that large values of $\epsilon$ help reduce the variance of the estimator but increase its bias, while small values behave exactly opposite. Furthermore, if we compare Eq. 9 to the mean squared error (MSE) of the empirical distribution $\widetilde{F}^{N}(t;x)$, which equals $F(t;x)\left[1-F(t;x)\right]/N$, we conclude that our estimation has a smaller MSE than the empirical distribution if $\epsilon$ is chosen optimally. Azzalini in [2] also shows that the properties of $\hat{F}^{N}_{\epsilon}$ extend to the quantile estimator, i.e., to the root of $\hat{F}^{N}_{\epsilon}(t;x)=1-\alpha$.

These observations regarding the behavior of the quantile approximation coincide with the findings of our numerical experiments. The solutions returned by the solver for large values of $\epsilon$ are generally feasible but sub-optimal due to the bias introduced by $\epsilon$. Large values of $\epsilon$ usually generate a conservative approximation of the quantile function; however, the solutions have less variance. We also find that solving the approximation with larger values of $\epsilon$ tends to require fewer iterations. Thus, there is a trade-off in the choice of $\epsilon$; small values of $\epsilon$ reduce bias, whereas larger values of $\epsilon$ reduce the variance and the computational effort. In section 6.1, we suggest a method for selecting an appropriate value of $\epsilon$ for an instance. Our proposal to find a positive $\epsilon$ that balances bias and variance contrasts with prior work [8, 15, 25, 26] that uses smoothing approximations which focused on driving $\epsilon$ to zero.

For the experiments in Section 6 we use the quartic kernel [39, p. 353] to define function $\Gamma_{\epsilon}$. This gives the following $\gamma_{\epsilon}$;

[TABLE]

which makes $\Gamma_{\epsilon}(\cdot)$ a twice continuously differentiable function.

Note 2.4**.**

*For simplicity of notation, we assume in the following sections that the quantile of interest is always the $(1-\alpha)$-quantile and drop the superscript of the quantile approximation. Therefore, the $(1-\alpha)$-quantile $\epsilon$-approximation is denoted just as $Q_{\epsilon}$ for now onward. *

3 Convergence and feasibility of the approximation

In this section we study how the solutions obtained by our method approximate those of the real problem Eq. 1. Throughout the rest of the paper we make the following assumption.

{assumption}

For each $x\in X$, the random variable $C(x,\xi)$ has a continuous distribution.

Notice that, if $C(x,\xi)$ is not a continuous random variable, we can make Section 3 hold for an approximation of the problem by adding a continuous random variable $\bar{\xi}$ to $C(x,\xi)$ with expectation zero and small variance. This could be treated as a fixed approximation, or one could solve a sequence of approximations in which the variance of $\bar{\xi}$ is decreased to zero.

Let the feasible region of the original problem Eq. 1 be defined as

[TABLE]

where $F(t;x)=\mathbb{P}(C(x,\xi)\leq t)$. Moreover, define $F_{\epsilon}(t;x)$, for $\epsilon>0$, as

[TABLE]

where $\Gamma_{\epsilon}$ is the smooth approximation of the indicator function given by Eq. 6.

Now, let $\{\xi_{i}\}_{i=1}^{N}$ be a random sample of the variable $\xi$ of size $N$. As in section 2.1, we approximate $F_{\epsilon}$ with $F_{\epsilon}^{N}$, and denote the feasible region defined by the approximation as follows

[TABLE]

Notice that if $\delta=\alpha$ and $t=0$, then $X_{\epsilon,\alpha}^{N,0}$ corresponds to the feasible region of Eq. 8. We consider this approximation more generally with $\delta\leq\alpha$ and $t\geq 0$ as we will show how the parameters $t$ and $\delta$ in Eq. 12 can be chosen to increase the likelihood that the solution to the approximate model is feasible to the true model. First, notice that if $t=\epsilon$, then $F_{\epsilon}^{N}(-t;x)\leq\widetilde{F}^{N}(0;x)$. This implies that a point $x\in X_{\epsilon,\alpha}^{N,\epsilon}$ satisfies $\widetilde{F}^{N}(0;x)\geq 1-\alpha$, so the solutions returned by our method are feasible for the empirical approximation of the feasible region. Moreover, the $\delta$ parameter can be used, as in [29], to ensure that the solutions obtained are feasible with high probability.

We start the analysis of Eq. 8 by showing asymptotic convergence of the optimal value of the approximation to the optimal value of the true problem. In particular, we analyze convergence of the approximation with respect to both the smoothing parameter and the sample size, whereas previous results, to the best of our knowledge, only consider convergence with respect to the smoothing parameter ([15, 25, 26]) or the sample size ([29, 36]) separately. Then, we show complementary analysis of the approximation problem that indicates that, with appropriate choice of parameters, any feasible solution to the approximation problem will be feasible with high probability for sufficiently large sample size.

3.1 Asymptotic convergence

In the following theorems we state results analogous to Proposition 2.1 and 2.2 in [36], thus the proofs we present are similar to those in that paper. We make the following assumptions.

{assumption}

$C(x,\xi)$ is a Carathéodory function, i.e., $C(x,\cdot)$ is measurable for every $x\in X$ and $C(\cdot,\xi)$ continuous for a.e. $\xi\in\Xi$.

{assumption}

There exists a globally optimal solution of Eq. 1, $\bar{x}$, such that for any $\delta>0$ there is $x\in X$ such that $\|x-\bar{x}\|\leq\delta$ and $F(0;x)>1-\alpha$.

This assumption implies that there exists a sequence $\{x_{k}\}_{k=1}^{\infty}\subseteq X$ that converges to an optimal solution $\bar{x}$ such that $F(0;x_{k})>1-\alpha$ for all $k\in\mathbb{N}$.

Theorem 3.1**.**

Suppose that Sections 3 and 3.1 hold. Then, $F_{\epsilon}^{N}(0;x)$ converges to $F(0;x)$ uniformly on any compact set $U\subseteq X$ w.p.1, i.e.,

[TABLE]

Proof 3.2**.**

We prove uniform convergence of $F_{\epsilon}^{N}(0;x)$ to $F(0;x)$ by first showing uniform convergence of $F_{\epsilon}^{N}(0;x)$ to $F_{\epsilon}(0;x)$ over $x\in U$ and $\epsilon\in[0,1]$ as $N\rightarrow\infty$, and then uniform convergence of $F_{\epsilon}(0;x)$ to $F(0;x)$ over $x\in U$ as $\epsilon\rightarrow 0$.

To prove uniform convergence of $F_{\epsilon}^{N}(0;x)$ to $F_{\epsilon}(0;x)$ in the compact set $U\subseteq X$ for all $\epsilon\in[0,1]$, we start by showing that $\Gamma_{\epsilon}(C(x,\xi))$ is continuous w.p.1 for all $x\in U$ and $\epsilon\in[0,1]$. Notice that $\Gamma_{\epsilon}(C(x,\xi))$ is continuous w.p.1 for all $x\in U$ and $\epsilon>0$. This follows because $\Gamma_{\epsilon}(\cdot)$ is continuous and $C(x,\xi)$ is a Carathéodory function. It remains to prove that $\Gamma_{\epsilon}(C(x,\xi))$ is continuous in $x$ and $\epsilon$ at $\epsilon=0$ w.p.1. Let $x$ and $\xi$ be such that $C(x,\xi)$ is continuous and strictly positive. Then, there exists $\epsilon_{0}>0$ and $\delta_{0}>0$ such that $C(x,\xi)>\epsilon_{0}$ and $C(y,\xi)>\epsilon_{0}$ for all $y\in\mathcal{B}_{\delta_{0}}(x)$ (where $\mathcal{B}_{\delta_{0}}(x)=\{y\mid\|x-y\|\leq\delta_{0}\}$). Thus, for all $\rho>0$, $|\Gamma_{0}(C(x,\xi))-\Gamma_{\epsilon}(C(y,\xi))|=0<\rho$ for all $0\leq\epsilon\leq\epsilon_{0}$ and $y\in\mathcal{B}_{\delta_{0}}(x)$. This shows that $\Gamma_{0}(C(x,\xi))$ is continuous when $C(x,\xi)>0$. An analogous proof can be derived if $C(x,\xi)$ is strictly negative. Finally, because $C(x,\xi)\neq 0$ w.p.1 for all $x\in X$ from Section 3, we conclude that $\Gamma_{\epsilon}(C(x,\xi))$ is continuous w.p.1 for all $x\in U$ and $\epsilon\in[0,1]$.

Uniform convergence of $F_{\epsilon}^{N}(0;x)$ to $F_{\epsilon}(0;x)$ over $x\in U$ and $\epsilon\in[0,1]$ is then a result of applying Theorem 7.53 in [40]. This is because function $\Gamma_{\epsilon}(C(x,\xi))$ is continuous w.p.1 for all $x\in U$ and $\epsilon\in[0,1]$, $|\Gamma_{\epsilon}(C(x,\xi))|\leq 1$ for all $x\in X$, and the sample is i.i.d.

Uniform convergence of $F_{\epsilon}(0;x)$ to $F(0;x)$ in the compact set $U\subseteq X$ as $\epsilon\rightarrow 0$ is a result of applying Theorem 4.91 in [40]. This comes from Sections 3 and 3.1.

Combining these two results, for every $\delta>0$, we can find $\epsilon_{\delta}>0$ and $N_{\delta}\in\mathbb{N}$ such that $|F_{\epsilon}(0;x)-F(0;x)|<\delta/2$ and $|F_{\epsilon}^{N}(0;x)-F_{\epsilon}(0;x)|<\delta/2$, w.p.1, for all $x\in U$, $\epsilon\leq\min\{1,\epsilon_{\delta}\}$, and $N\geq N_{\delta}$, which ensures that

[TABLE]

*We can therefore conclude that $F_{\epsilon}^{N}(0;x)$ converges to $F(0;x)$ uniformly in the compact set $U$. *

We denote the sets of optimal solutions for problems Eq. 1 and Eq. 8 as $S$ and $S_{\epsilon}^{N}$, respectively. We also denote $v^{*}$ and $v_{\epsilon}^{N}$ the optimal values of Eq. 1 and Eq. 8, respectively. We show next that, under some assumptions, $v^{N}_{\epsilon}$ and $S^{N}_{\epsilon}$ converge w.p.1 to their counterparts of the true problem as $N$ increases and $\epsilon$ decreases.

Theorem 3.3**.**

*Suppose that $X$ is compact, the function $f$ is continuous, and Sections 3, 3.1, and 3.1 hold. Then, $v_{\epsilon}^{N}\rightarrow v^{*}$ and $\mathbb{D}(S_{\epsilon}^{N},S)\rightarrow 0$ w.p.1 as $N\rightarrow\infty$ and $\epsilon\rightarrow 0$. *

Proof 3.4**.**

Notice that for any given $N$ and $\epsilon$, a point $x$ in the feasible set of Eq. 8 satisfies $F_{\epsilon}^{N}(0;x)\geq 1-\alpha$. Thus, we can show convergence of the feasible set of Eq. 8 to the feasible set of Eq. 1 by analyzing the function $F_{\epsilon}^{N}$.

Now, by section 3.1, the set $S$ is nonempty and there is $x\in X$ such that $F(0;x)>1-\alpha$. Because $F_{\epsilon}^{N}(0;x)$ converges to $F(0;x)$, by Theorem 3.1, there exists $\epsilon_{0}$ small enough and $N_{0}$ large enough such that $F_{\epsilon_{0}}^{N_{0}}(0;x)\geq 1-\alpha$ w.p.1. Because $F_{\epsilon_{0}}^{N_{0}}$ is continuous in $x$ and $X$ is compact, the feasible set of the approximation problem Eq. 8 is compact as well, and hence $S_{\epsilon}^{N}$ is nonempty w.p.1 for all $N\geq N_{0}$ and $\epsilon\leq\epsilon_{0}$.

Let $\{N_{k}\}_{k=1}^{\infty}\geq N_{0}$ and $\{\epsilon_{k}\}_{k=1}^{\infty}\leq\epsilon_{0}$ be two sequences such that $N_{k}\rightarrow\infty$ and $\epsilon_{k}\rightarrow 0$. Let $\hat{x}_{k}\in S_{\epsilon_{k}}^{N_{k}}$, i.e., $\hat{x}_{k}\in X$, $F_{\epsilon_{k}}^{N_{k}}(0;\hat{x}_{k})\geq 1-\alpha$ and $v_{\epsilon_{k}}^{N_{k}}=f(\hat{x}_{k})$. Let $\hat{x}\in X$ be any cluster point of $\{x_{k}\}_{k=1}^{\infty}$ and $\{x_{l}\}_{l=1}^{\infty}$ be a subsequence converging to $\hat{x}$. Because $F_{\epsilon_{l}}^{N_{l}}(0;x)$ is continuous and converges uniformly to $F(0;x)$ on $X$ w.p.1 by Theorem 3.1, we have that $F(0;\hat{x})=\lim_{l\rightarrow\infty}F_{\epsilon_{l}}^{N_{l}}(0;\hat{x}_{l})$, w.p.1. Hence, $F(0;\hat{x})\geq 1-\alpha$ and $\hat{x}$ is feasible for the true problem, and therefore $f(\hat{x})\geq v^{*}$. Also $f(\hat{x}_{l})\rightarrow f(\hat{x})$ w.p.1 which means that $\lim_{l\rightarrow\infty}v_{\epsilon_{l}}^{N_{l}}\geq v^{*}$ w.p.1. Since this is true for any cluster point of $\{x_{k}\}_{k=1}^{\infty}$ in the compact set $X$, we have

[TABLE]

Now, by section 3.1, there exists an optimal solution $\bar{x}$ and a sequence $\{x_{l}\}_{l=1}^{\infty}$ converging to $\bar{x}$ with $F(0;x_{l})>1-\alpha$. Since $F_{\epsilon_{k}}^{N_{k}}(0;{x}_{l})$ converges to $F(0;x_{l})$ w.p.1, there exist $K(l)$ such that $F_{\epsilon_{k}}^{N_{k}}(0;{x}_{l})\geq 1-\alpha$ for every $k\geq K(l)$ and every $l$, w.p.1. Without loss of generality we can assume that $K(l)<K(l+1)$ for every $l$ and define the sequence $\{\bar{x}_{k}\}_{k=K(1)}^{\infty}$ by setting $\bar{x}_{k}=x_{l}$ for all $k$ and $l$ with $K(l)\leq k<K(l+1)$. We then have $F_{\epsilon_{k}}^{N_{k}}(0;\bar{x}_{k})\geq 1-\alpha$, which implies $v_{\epsilon_{k}}^{N_{k}}\leq f(\bar{x}_{k})$ for all $k\geq K(1)$. Since $f$ is continuous, we have that

[TABLE]

*It follows from Eq. 13 and Eq. 14 that $v_{\epsilon_{k}}^{N_{k}}\rightarrow v^{*}$ w.p.1. To show that $\mathbb{D}(S_{\epsilon}^{N},S)\rightarrow 0$ w.p.1. we defer to Theorem 5.3 in [40] which can easily be adapted to our context. *

Notice that Theorem 3.3 only applies to the globally optimal objective function and solutions.

3.2 Finite sample feasibility of approximation problem solutions

We now consider conditions under which an optimal solution for the approximation problem Eq. 8, if it exists, is feasible for the true problem Eq. 1. These results generalize those presented in Section 2.2 in [29], which considers the empirical cdf instead of a smooth approximation; therefore, some of our proofs follow from the ones presented there.

We will make use of Hoeffding’s inequality.

Theorem 3.5** (Hoeffding’s Inequality [24]).**

Let $Y_{1},\ldots,Y_{N}$ be independent random variables, with $\operatorname{\mathbb{P}}(Y_{i}\in[a_{i},b_{i}])=1$, where $a_{i}\leq b_{i}$ for $i=1,\ldots,N$. Then if $t>0$,

[TABLE]

{assumption}

Let the scalars $t\in\mathbb{R}$ and $\delta\in[0,\alpha]$ be chosen such that $M:=\inf_{x\in X}M_{x}>0$, where $M_{x}=F(0;x)-F_{\epsilon}(-t;x)+(\alpha-\delta)$.

The above assumption can always be satisfied if we pick $t=\epsilon$ and $\delta<\alpha$. This claim can easily be observed from the definition of $F_{\epsilon}(-t;x)$ in Eq. 11, which implies that $F_{\epsilon}(-\epsilon;x)\leq F(0;x)$ for all $x\in X$.

Theorem 3.6** (Probabilistic feasibility guarantee).**

Let $x\in X$ be such that $x\notin X_{\alpha}$, and suppose that section 3.2 holds. Then

[TABLE]

Proof 3.7**.**

Let $Y_{i}=\Gamma_{\epsilon}(C(x,\xi_{i})+t)$ for all $i=1,\ldots,N$, then $\mathbb{P}(Y_{i}\in[0,1])=1$ and $\operatorname{E}[Y_{i}]=F_{\epsilon}(-t;x)$. Given that $x$ is not a feasible solution for the original problem, $F(0;x)<1-\alpha$. We therefore have

[TABLE]

*where the last inequality is obtained from Hoeffding’s Inequality. *

Theorem 3.6 shows that if we shift the approximation of the cdf appropriately and/or decrease the risk level from $\alpha$ to $\delta\leq\alpha$, then a point $x$ that is feasible for the original problem Eq. 1 is also in $X_{\epsilon,\delta}^{N,t}$ with a probability that increases exponentially with the size of the sample $N$. Furthermore, under certain conditions, we can set $t=0$ and $\delta=\alpha$ and still get the results stated above. This can be seen in the following corollary.

Corollary 3.8**.**

Let $x\in X$ be such that $Y_{x}=C(x,\xi)$ is a continuous random variable with a strictly decreasing probability density function (pdf) $h_{x}(y)$ in the interval $(-\epsilon,\epsilon)$. If $F(0;x)<1-\alpha$, then

[TABLE]

where $\beta_{x}=F(0;x)-F_{\epsilon}(0;x)>0$.

Proof 3.9**.**

This follows from Theorem 3.6 for $t=0$, $\delta=\alpha$, and noticing that $\beta_{x}=M_{x}=F(0;x)-F_{\epsilon}(0;x)>0$ because

[TABLE]

*The fourth equality is the result of defining $\Gamma_{\epsilon}$ as the integral of a symmetric pdf, which implies that $1-\Gamma_{\epsilon}(y)=\Gamma_{\epsilon}(-y)$. The last inequality is due to the fact that the pdf of $Y_{x}$ is strictly decreasing in the interval $(-\epsilon,\epsilon)$, which implies that $h_{x}(-y)>h_{x}(y)$ for $y\in(0,\epsilon)$, and due to $\Gamma_{\epsilon}(y)=\gamma_{\epsilon}(y)>0$ for $y\in[0,\epsilon)$. *

One of the properties that follows from Corollary 3.8 is that, for $\epsilon>0$, our approximation is asymptotically conservative. In other words, for a fixed $\epsilon>0$ and under a monotonicity assumption, a point $x$ that is feasible for the original problem Eq. 1 is also feasible for the approximation problem Eq. 8 with a probability that increases exponentially with the sample size $N$. It can be shown that for the vast majority of the common continuous distributions such as exponential, normal and Weibull, to mention a few, there exists a neighborhood of the solution such that the assumption in Corollary 3.8 holds.

Corollary 3.10**.**

Assume that $\Gamma_{\epsilon}$ is defined as in Eq. 6, where $\gamma_{\epsilon}(\cdot)$ is given by Eq. 10. Let $x\in X$ be such that $Y_{x}=C(x,\xi)$ is a continuous random variable with pdf $h_{x}(y)$. If there exists $m_{1},m_{2}>0$ such that $h_{x}(y)\geq-m_{1}y+h_{x}(0)$ for all $y\in(-\epsilon,0]$, and $h_{x}(y)\leq-m_{2}y+h_{x}(0)$ for all $y\in[0,\epsilon)$. Then, we have that

[TABLE]

and thus, for all $x\notin X_{\alpha}$,

[TABLE]

Proof 3.11**.**

From the definition of $\beta_{x}=F(0;x)-F_{\epsilon}(0;x)$, we get

[TABLE]

*where $\int_{0}^{\epsilon}\Gamma_{\epsilon}(y)y\,dy=\int_{0}^{\epsilon}\gamma_{\epsilon}(y)y\,dy=\epsilon^{2}/28$. *

We now consider conditions under which the entire feasible region of the approximation problem is a subset of the feasible region of problem Eq. 1. For these results to hold we need to make additional assumptions. We begin with a preliminary result that holds under the assumption that the feasible region $X$ is finite.

Theorem 3.12**.**

Suppose that $X$ is finite, and that section 3.2 holds. Then

[TABLE]

Proof 3.13**.**

Using the result from Theorem 3.6, we get the following

[TABLE]

We now extend this result to the case in which $X$ is not finite. For this result, we make the following Lipschitz continuity assumption on $C$.

{assumption}

There exists $L>0$ such that

[TABLE]

This assumption was made in [29]. The proof stated below follows the same ideas as Theorem 10 in [29].

Theorem 3.14**.**

Assume that $X$ is bounded with diameter $D$ and section 3.2 holds. Also, suppose that section 3.2 holds with $t>0$ and $\delta>0$. Then, for any $\beta\in\left(0,\delta\right]$, the following is true:

[TABLE]

Proof 3.15**.**

Let $J=\lceil 1/\beta\rceil$. For $j=1,\ldots,J-1$, define

[TABLE]

and let $X_{J}=\{x\in X\mid(J-1)/J\leq F(0;x)\leq 1\}$. Following the arguments used to prove Theorem 10 in [29], we have that for each $j$ there exists a finite set $Z_{j}\subseteq X_{j}$ such that $|Z_{j}|\leq\lceil 2LD/t\rceil^{n}$ and that for all $x\in X_{j}$ there exists $z\in Z_{j}$ such that $\|x-z\|_{\infty}\leq t/L$.

Now define $Z=\cup_{j=1}^{J}Z_{j}$, and notice that $|Z|\leq\lceil 1/\beta\rceil\lceil 2LD/t\rceil^{n}$. Next, define $Z_{\alpha-\beta}=\{x\in Z\mid F(0;x)\geq 1-\alpha+\beta\}$ and

[TABLE]

Since $Z$ is finite, we can apply Theorem 3.12 to obtain

[TABLE]

*Consider now any $x\in X_{\epsilon,(\delta-\beta)}^{N,2t}$. Let $j\in\{1,\ldots,J\}$ be such that $x\in X_{j}$. By definition of $Z_{j}$, there exists $z\in Z_{j}$ such that $\|x-z\|_{\infty}\leq t/L$. By section 3.2, we know that $|C(x,\xi_{i})-C(z,\xi_{i})|\leq t$ w.p.1. Therefore, $C(z,\xi_{i})+t\leq C(x,\xi_{i})+2t$, which implies that $\Gamma_{\epsilon}(C(z,\xi_{i})+t)\geq\Gamma_{\epsilon}(C(x,\xi_{i})+2t)$ w.p.1, and thus $F_{\epsilon}^{N}(-t;z)\geq F_{\epsilon}^{N}(-2t;x)\geq 1-\delta+\beta$ w.p.1. We conclude that if $x\in X_{\epsilon,(\delta-\beta)}^{N,2t}$, then there exists $z\in Z_{\delta-\beta}^{N,t}$ w.p.1. Furthermore, by the definition of $X_{j}$ and because $Z_{j}\subseteq X_{j}$, we have $|F(0;x)-F(0;z)|\leq\beta$. If $Z_{\delta-\beta}^{N,t}\subseteq Z_{\alpha-\beta}$ then $F(0;z)\geq 1-\alpha+\beta$, which implies that $F(0;x)\geq F(0;z)-\beta\geq 1-\alpha$, and the result follows. *

4 Single Chance-Constrained Problems

A single chance constrained problem is one in which $m=1$, which means that $C(x,\xi)=\max_{j=1,\ldots,m}\{c_{j}(x,\xi)\}=c(x,\xi)$. Defining $C^{N}(x)=c^{N}(x)=[c(x,\xi_{1}),\ldots,c(x,\xi_{N})]^{T}$, problem Eq. 8 then becomes

[TABLE]

The advantage of the single-chance constrained case is that the approximated probabilistic constraint can be formulated in such a way that it results in a single smooth constraint (assuming $c$ is smooth in $x$) that can be handled by standard nonlinear optimization solvers and can be added to any nonlinear optimization problem. Furthermore, for a given $x$ and $\epsilon>0$, we can use the implicit function theorem to obtain the gradient of $q$ (and provide it to the solver) as follows

[TABLE]

which is well defined as shown in Lemma 2.1.

In fact, the following proposition shows that if we define $\gamma_{\epsilon}$ as Eq. 10, then function $Q_{\epsilon}(\cdot)$ has Lipschitz continuous gradients. Using the chain rule, this property can be extended to $q(x)$ as long as $c(x,\xi)$ has Lipschitz continuous gradients for all instances $\xi_{i}$.

Proposition 4.1**.**

*If $\gamma_{\epsilon}$ is defined as Eq. 10 and $(1-\alpha)N\notin\mathbb{Z}$, then the function $Q_{\epsilon}:\mathbb{R}^{N}\rightarrow\mathbb{R}$ has bounded Lipschitz continuous gradients. *

Proof 4.2**.**

Let $z\in\mathbb{R}^{N}$. It follows from Lemma 2.1 and the implicit function theorem that

[TABLE]

We begin by showing that $\sum_{j=1}^{N}\Gamma_{\epsilon}^{\prime}(z_{j}-Q_{\epsilon}(z))$ is bounded away from zero. We prove this claim by contradiction. Assume that there exists $z\in\mathbb{R}^{N}$ such that $|z_{j}-Q_{\epsilon}(z)|\geq\epsilon$ for all $j\in\{1,\ldots,N\}$. Then, from Eq. 6, $\Gamma_{\epsilon}(z_{j}-Q_{\epsilon}(z))\in\{0,1\}$ for all $j$. This implies that $\sum_{j=1}^{N}\Gamma_{\epsilon}(z_{j}-Q_{\epsilon}(z))\in\mathbb{N}$, which means that $\sum_{j=1}^{N}\Gamma_{\epsilon}(z_{j}-Q_{\epsilon}(z))\neq(1-\alpha)N$. Therefore, there must exists $i\in\{1,\ldots,N\}$ and $0<\epsilon_{0}<\epsilon$ such that $|z_{i}-Q_{\epsilon}(z)|\leq\epsilon_{0}$ and $\sum_{j=1}^{N}\Gamma_{\epsilon}^{\prime}(z_{j}-Q_{\epsilon}(z))\leq\Gamma_{\epsilon}^{\prime}(z_{i}-Q_{\epsilon}(z))\leq\gamma_{\epsilon}^{\prime}(\epsilon_{0})<0$, for all $z\in\mathbb{R}^{N}$.

*More over, since $\sum_{j=1}^{N}\Gamma_{\epsilon}^{\prime}(z_{j}-Q_{\epsilon}(z))\leq\Gamma^{\prime}_{\epsilon}(z_{i}-Q_{\epsilon}(z))\leq 0$ for all $i$, we conclude that $0\leq[\nabla Q_{\epsilon}(z)]_{i}\leq 1$ for all $i$; this means that $\nabla Q_{\epsilon}(z)$ is bounded. Finally, because $\Gamma_{\epsilon}^{\prime\prime}$ exists and is bounded, then $\Gamma_{\epsilon}^{\prime}$ is Lipschitz continuous. Given that the denominator is bounded above by $\gamma_{\epsilon}^{\prime}(\epsilon_{0})<0$, we can conclude that $Q_{\epsilon}(z)$ has Lipschitz continuous gradients. *

If we also assume that $f$ and $c$ are twice continuously differentiable, then we can provide the Hessian of constraint Eq. 16b to the solver as well. The Hessian is then computed as

[TABLE]

where $\left[\nabla Q_{\epsilon}(\cdot)\right]_{i}$ is the $i$th entry of the gradient of $Q_{\epsilon}(\cdot)$, and $\nabla^{2}Q_{\epsilon}\left(\cdot\right)$ is the Hessian of the smooth approximation of the quantile with respect to each scenario, i.e., $\nabla^{2}Q_{\epsilon}\left(\cdot\right)\in\mathbb{R}^{N\times N}$ is such that

[TABLE]

where $w_{\ell}=\Gamma_{\epsilon}^{\prime}(z_{\ell}-Q(z))$, $w_{\ell}^{\prime}=\Gamma_{\epsilon}^{\prime\prime}(z_{\ell}-Q(z))$, $W=\sum_{\ell=1}^{N}w_{\ell}$, and $W^{\prime}=\sum_{\ell=1}^{N}w_{\ell}^{\prime}$.

In section 6.2 we discuss the performance of this approximation. We give two examples in which the nonlinear constraint $q(x)$ and its derivatives are given directly to Knitro.

5 Joint Chance-Constrained Problems

In this section, we consider problems where $c(x,\xi)$ is a vector-valued function from $\mathbb{R}^{n}\times\mathbb{R}^{s}$ to $\mathbb{R}^{m}$ for $m>1$, and therefore $C(x,\xi)=\max_{j=1,\ldots,m}\{c_{j}(x,\xi)\}$ is not smooth even when $c_{j}(x,\xi)$ are smooth for $j=1,\ldots,m$.

Assuming that there exists a deterministic differentiable function, $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$, such that the deterministic feasible set $X$ can be written as $X=\left\{x\in\mathbb{R}^{n}\mid g(x)\leq 0\right\}$, we can rewrite the approximation of Eq. 8 as

[TABLE]

The discussion in this section can be extended to solve problems with deterministic equality constraints in the definition of the set $X$, but we consider only inequality constraints for a succinct exposition.

Notice that, because $C(x,\xi)$ is no longer smooth, the constraint (19c) is not necessarily differentiable. To obtain a smooth reformulation of the quantile constraint Eq. 19c, we introduce auxiliary variables $z_{i}$ that correspond to $\max_{j}\{c_{j}(x,\xi_{i})\}$ and consider

[TABLE]

where $\boldsymbol{1}_{m}$ is a vector of ones of size $m$.

Proposition 5.1**.**

*Let $x^{*}$ be a local minimum of Eq. 19, then $(x^{*},C^{N}(x^{*}))$ is a local minimum of Eq. 20. *

Proof 5.2**.**

*Assume that $(x^{*},C^{N}(x^{*}))$ is not a local minimum of Eq. 20. This means that there exists $(\hat{x},\hat{z})$ in a neighborhood of $(x^{*},C^{N}(x^{*}))$ such that $f(\hat{x})<f(x^{*})$, $g(\hat{x})\leq 0$, $c_{j}(\hat{x},\xi_{i})\leq\hat{z}_{i}$ for all $i$ and $j$, and $Q_{\epsilon}(\hat{z})\leq 0$; thus, $[C(\hat{x},\xi)]_{i}\leq\hat{z}_{i}$ for all $i$. Now, the function $Q_{\epsilon}$ is an increasing function with respect to $\hat{z}_{i}$ because $0\leq[\nabla Q_{\epsilon}(\hat{z})]_{i}$, as shown in Proposition 4.1. This means that $Q(C^{N}(\hat{x},\xi))\leq Q(\hat{z})\leq 0$. Putting everything together, we see that $\hat{x}$ is feasible for (5.1) and $f(\hat{x})<f(x^{*})$. This argument holds for $(\hat{x},\hat{z})$ from an arbitrarily small neighborhood of $(x^{*},C^{N}(x^{*}))$. This is a contradiction because $x^{*}$ is assumed to be a local minimum of (5.1). *

We have shown that a local minimum of Eq. 19 is also a local minimum of Eq. 20 if we define $z_{i}^{*}=[C^{N}(x^{*})]_{i}$. However, a local minimum of Eq. 20 is not necessarily a local minimum of Eq. 19 if $z_{i}^{*}\neq[C^{N}(x^{*})]_{i}$. Thus, the downside of formulating the JCCP as the SCCP in Eq. 20 is that a standard NLP solver might get stuck at a spurious local minimum. We therefore introduce a tailored algorithm that is designed to solve JCCPs.

5.1 Exact-penalty formulation

In order to solve problem Eq. 19, we write an equivalent unconstrained optimization problem in which the constraints are handled by an $\ell_{1}$-penalty term. Our approach is motivated by the S$\ell_{1}$QP method for regular NLPs [35]. We define the penalty function

[TABLE]

where $\pi>0$ is a penalty parameter and $[x]^{+}_{i}=\max\{0,x_{i}\}$. If $x^{*}$ is a minimizer of $\phi_{\pi}$ for $\pi>0$ and $\|[g(x)]^{+}\|_{1}+[Q_{\epsilon}(C^{N}(x))]^{+}=0$, then $x^{*}$ solves Eq. 19 [23, p. 299]. The following theorem, stated in [11, Theorem 2.1], considers the converse situation and shows that there exists $\pi>0$ such that we can use the minimization of (21) to obtain a solution for Eq. 19.

Theorem 5.3**.**

Suppose $f$, $g$, $c$, and $Q_{\epsilon}$ are locally Lipschitz on $\mathbb{R}^{n}$ and let $x^{*}$ be a local minimizer of Eq. 19. Moreover, suppose that the Mangasarian-Fromowitz condition holds for problem Eq. 20 at $x^{*}$ and $z^{*}=C^{N}(x^{*})$. Then, there exists $\pi^{*}>0$ and $\delta>0$ such that for all $x\in\mathbb{B}_{\delta}(x^{*})$,

[TABLE]

To minimize function Eq. 21, we propose a trust-region method that generates steps from the minimization of a piecewise quadratic model of $\phi_{\pi}$ in every iteration. By making sure that the model approximates $\phi_{\pi}$ to first-order, we can invoke existing convergence results that show that the trust-region method converges to a stationary point of $\phi_{\pi}$. A piecewise quadratic model with such characteristics is defined in the following proposition.

Proposition 5.4**.**

Let $f$ and $g$ be differentiable with Lipschitz continuous gradients. Also, given a sample $\{\xi_{1},\ldots,\xi_{N}\}$, let $c(x,\xi_{i})$ be differentiable with respect to $x$ for all $i\in\{1,\ldots,N\}$ and assume there exists $L>0$ and $L^{\prime}>0$ such that

[TABLE]

Given a point $x\in\mathbb{R}^{n}$ and a direction $d\in\mathbb{R}^{n}$, consider the function

[TABLE]

where $H\in\mathbb{R}^{n\times n}$ is a symmetric matrix, and

[TABLE]

Then, $m(x,H;d)$ is a first-order model of $\phi_{\pi}(x+d)$ in the sense that

[TABLE]

for all $d,x\in\mathbb{R}^{n}$ and all $H$ with $\|H\|\leq L_{H}$, for fixed $L_{H}$.

Proof 5.5**.**

We only prove that $|Q_{\epsilon}(C^{N}(x+d))-\widetilde{Q}_{\epsilon,C}(x;d)|=\mathcal{O}(\|d\|^{2})$ since the rest follows from the Lipschitz continuity of $\nabla f$ and $\nabla g$, and boundedness of $H$.

For brevity in the proof, we will use the following notation $C_{y}=C^{N}(y)$ and $\widetilde{C}_{(x;d)}=\widetilde{C}^{N}(x;d)$. We have that,

[TABLE]

Analyzing the first term, from Taylor’s Theorem we know that there exists $\zeta_{1}\in(0,1)$ such that

[TABLE]

where the two last inequalities come from Proposition 4.1 and Eq. 22 (we denote $L_{Q}$ as the Lipschitz constant of $\nabla Q_{\epsilon}$).

For the second term we have

[TABLE]

Because $\nabla Q_{\epsilon}$ by Proposition 4.1 is bounded, it only remains to show that $\|C_{x+d}-\widetilde{C}_{(x;d)}\|=\mathcal{O}(\|d\|^{2})$.

To prove this, first notice that

[TABLE]

Therefore, from Taylor’s Theorem we have that there exists $\zeta_{2}\in(0,1)$ such that

[TABLE]

where $j_{1}^{i}\in\operatorname{argmax}_{j}\{c_{j}(x+d,\xi_{i})\}$. Analogous to Eq. 25, we get

[TABLE]

*where $j_{2}^{i}\in\operatorname{argmax}_{j}\{c_{j}(x,\xi_{i})+\nabla c_{j}(x,\xi_{i})^{T}d\}$. It follows from Eq. 25 and Eq. 26 that $\|C_{x+d}-\widetilde{C}_{(x;d)}\|\leq\sqrt{N}L^{\prime}\|d\|^{2}$. *

Now that model Proposition 5.4 has been established, we proceed to give an outline of the algorithm used to minimize $\phi_{\pi}$.

5.2 Trust-Region Algorithm

In this section we present an S$\ell_{1}$QP-type trust-region algorithm to minimize Eq. 21. This algorithm seeks to obtain a first-order stationary point of $\phi_{\pi}$ by solving a sequence of quadratic programs (QP). At each iteration $k$, we compute a step $d^{k}$ by minimizing the model $m(x^{k},H^{k};d)$ within a radius $\Delta_{k}$ for a given $H_{k}$. This problem is equivalent to

[TABLE]

We describe the choice of $H^{k}$ in section 5.3.

A step $d^{k}$, obtained from Eq. 27, is accepted if it results in sufficient decrease of $\phi_{\pi}$, i.e., the step $d^{k}$ is taken only if the value $\phi_{\pi}(x^{k}+d^{k})$ is sufficiently smaller than $\phi_{\pi}(x^{k})$. We then proceed with iteration $k+1$, where now $x^{k+1}=x^{k}+d^{k}$. Algorithm 1 states the full algorithm.

The following theorem follows from [10, Theorem 11.2.5]. It states that any limit point $x^{*}$ of the sequence $\{x^{k}\}$ generated by Algorithm 1 is a stationary point of function $\phi_{\pi}$. Thus, we conclude that $x^{*}$ is a local minimum of Eq. 19 if $x^{*}$ is a local minimizer of $\phi_{\pi}$, $g(x^{*})\leq 0$, and $Q_{\epsilon}\left(C^{N}(x^{*})\right)\leq 0$.

Theorem 5.6**.**

*Suppose that Algorithm 1 does not terminate at a stationary point of $\phi_{\pi}$ in 3, and that $x^{*}$ is a limit point of the sequence $\{x^{k}\}$ generated by Algorithm 1. Let the assumptions from Proposition 5.4 hold and assume that, at each iteration, problem Eq. 27 is solved exactly. Then $x^{*}$ is a first-order critical point of $\phi_{\pi}$. *

To apply [10, Theorem 11.2.5], it is necessary that $m(x,H;d)$ is a first-order model, which was shown in Proposition 5.4. The condition that problem Eq. 27 is solved exactly can be relaxed as long as some Cauchy decrease condition holds (refer to [10, Theorem 11.2.5]).

5.3 Choice of Hessian

In this section we propose a choice of the Hessian $H^{k}$ in Eq. 27 that aims to attain fast local convergence. The key observation is a close relationship between the non-smooth formulation Eq. 19 and the smooth optimization problem

[TABLE]

where $\overline{C}^{N}(x)$ is defined below in Lemma 5.7.

Recall that, if $x^{*}$ is a local minimizer of Eq. 19, then $x^{*}$ and $z^{*}=C^{N}(x^{*})$ is a KKT point of Eq. 20 under some constraint qualification (see Theorem 5.3). The following lemma shows that $x^{*}$ is also a KKT point of Eq. 28.

Lemma 5.7**.**

Let $(x^{*},z^{*},\nu^{*},\mu^{*},\lambda^{*})$ be a KKT point of Eq. 20 with $z^{*}=C^{N}(x^{*})$, where $\nu^{*}$, $\mu^{*}$ and $\lambda^{*}$ are the multipliers associated to constraints Eq. 20b, Eq. 20c and Eq. 20d, respectively. For each $i\in\{1,\ldots,N\}$, set

[TABLE]

if $\lambda^{*}\neq 0$ and $\left[\nabla Q_{\epsilon}(C^{N}(x^{*}))\right]_{i}\neq 0$. Else, select one $j$ such that $c_{j}(x^{*},\xi_{i})=C(x^{*},\xi_{i})$ and define $[\bar{\mu}_{i}]_{j}=1$ and $[\bar{\mu}_{i}]_{\ell}=0$ if $\ell\neq j$. Also, define the following function

[TABLE]

*Then, $(x^{*},\nu^{*},\lambda^{*})$ is a KKT point of Eq. 28. *

Proof 5.8**.**

Because $(x^{*},z^{*},\nu^{*},\mu^{*},\lambda^{*})$ is a KKT point of Eq. 20, the following conditions hold

[TABLE]

where $\nabla g$ and $\nabla c$ correspond to the transpose of the Jacobians of $g$ and $c$, respectively. Notice that the KKT conditions of Eq. 20 associated with $g$ have been left out because they trivially extend to the KKT conditions of Eq. 28.

If $\lambda^{*}\neq 0$ and $\left[\nabla Q_{\epsilon}(C^{N}(x^{*}))\right]_{i}\neq 0$, then

[TABLE]

If $\lambda^{*}=0$ or $\left[\nabla Q_{\epsilon}(C^{N}(x^{*}))\right]_{i}=0$, it follows from the definition of $\bar{\mu}_{i}$, that $[\overline{C}^{N}(x^{*})]_{i}=\bar{\mu}_{i}^{T}c(x^{*},\xi_{i})=z_{i}^{*}$. Thus, feasibility, $Q_{\epsilon}(\overline{C}^{N}(x^{*}))=Q_{\epsilon}(z^{*})\leq 0$, and complementarity, $\lambda^{*}Q_{\epsilon}(z^{*})=\lambda^{*}Q_{\epsilon}(\overline{C}^{N}(x^{*}))=0$, hold.

Now, note that the gradient of the Lagrangian of Eq. 28 is given by

[TABLE]

Then,

[TABLE]

*Thus, $(x^{*},\nu^{*},\lambda^{*})$ is a KKT point of Eq. 28. *

Lemma 5.7* shows that we can use the norm of Eq. 33 in the termination criterion for our algorithm. Also, from complementary slackness and feasibility, we have that most of the entries of the vector $\mu_{i}^{*}$ are zero. The reason is that only the $j$th entries of $\mu_{i}^{*}$ such that $c_{j}(x^{*},\xi_{i})=C(x^{*},\xi_{i})$ can be nonzero. This means that the positive entries of $\mu_{i}^{*}$ indicate which entries of vector $c(x^{*},\xi_{i})$ correspond to the maximum $C(x^{*},\xi_{i})$. Therefore, $\overline{C}^{N}(\cdot)$ can be interpreted as a smooth function that coincides with the maximum $C^{N}(\cdot)$ at the solution $x^{*}$.*

Suppose we apply the basic SQP algorithm [35] to Eq. 28 and compute steps $d^{k}$ from the solution of the QP

[TABLE]

to update the iterate $x^{k+1}=x^{k}+d^{k}$. If ${H}^{k}$ is chosen as the Hessian of the Lagrangian function $\overline{\mathcal{L}}$ of Eq. 28 at the current iterate $(x^{k},\nu^{k},\lambda^{k})$, i.e.,

[TABLE]

where $[\nabla\overline{C}^{N}(x^{k})]_{\boldsymbol{\cdot}i}=\nabla c(x,\xi_{i})\bar{\mu}_{i}$ for all $i\in\{1,\ldots,N\}$, then this algorithm exhibits quadratic local convergence, under standard regularity assumptions.

Our goal is to mimic this behavior in Algorithm 1. Using arguments similar to those in the proof of Lemma 5.7 one can show that the steps generated by the QP Eq. 27 converge to the SQP steps from Eq. 34 up to first order as $x^{k}\to x^{*}$, under suitable regularity assumptions. This suggests to choose $H^{k}$ defined by Eq. 35 for the step computation Eq. 27 in Algorithm 1, with the hope that the steps used in Algorithm 1 approach the SQP steps for Eq. 28 and result in fast local convergence.

However, the definition of Eq. 28 requires the knowledge of the optimal solution in Eq. 29, which is not available. Instead, we approximate the optimal solution with the current iterate and define

[TABLE]

Then, the Hessian matrix in the QP Eq. 27 is chosen as follows:

[TABLE]

where $[\nabla\overline{C}^{N}(x^{k})]_{\boldsymbol{\cdot}i}=\nabla c(x,\xi_{i})\bar{\mu}^{k}_{i}$ for all $i\in\{1,\ldots,N\}$. Note the subtle difference with (35): The argument at which derivatives of $Q_{\epsilon}$ are evaluated is based on the true maximum values given by the non-smooth function $C^{N}(x^{k})$ instead of the smoothed version $\overline{C}^{N}(x^{k})$. With this choice of $H^{k}$, we observed superlinear local convergence of Algorithm 1 in our experiments.

6 Numerical results

In this section we report results demonstrating the performance of our proposed approach on two single chance-constrained (SCC) examples and one joint chance-constrained (JCC) example. The goal is to observe the performance of our method from different perspectives such as: solution quality, speed, and features related to the smoothing parameter $\epsilon$.

In all the experiments we use the $\Gamma_{\epsilon}$ function defined from $\gamma_{\epsilon}$ in Eq. 10, which defines a twice continuously differentiable function. This allows us to incorporate first- and second-order derivatives into the algorithms. All computations were executed on Ubuntu 16.04 with 256GB RAM and a CPU with two Intel Xeon processors with 10 cores each, running at 3.10GHz. The algorithm and all experiments were implemented in Matlab R2015b.

In these experiments, because our choices of $\alpha$ and $N$ are such that $(1-\alpha)N\in\mathbb{Z}$, we opt to add a constant $b=1/2$ to the left-hand side of Eq. 7 in our implementation (refer to Remark 2.3).

6.1 Choosing the smoothing parameter

We introduced an approximation of problem Eq. 1 that does not need explicit information about the distribution of the random variables, their mean or their variance, and that only requires a sample to compute the constraint approximation. However, the quality of the solutions obtained from the approximation depends on the choice of the parameter $\epsilon$. As established in Section 3, $\epsilon$ is related to the feasibility of the solutions. Moreover, in our experiments we observed that large values of $\epsilon$ are associated with fewer iterations and less variability in the solutions but at the same time lead to conservativeness, while small values work in the opposite manner. Therefore, a good choice of $\epsilon$ is of vital importance to the solution quality.

For a fixed sample, we propose a binary search algorithm to find a value of $\epsilon$ which yields a solution that is feasible to the chance constraint, but not too conservative. We begin by setting $\ell=0$, selecting an initial $\epsilon_{0}$, and setting $\epsilon_{LB}=0$ and $\epsilon_{UB}=\infty$. Then, we solve problem Eq. 19 for $\epsilon_{\ell}$ and estimate the true probability $p_{\ell}=\operatorname{\mathbb{P}}(C(x_{\ell}^{*},\xi)\leq 0)$, where $x_{\ell}^{*}$ is the optimal solution of problem Eq. 19 using the parameter $\epsilon_{\ell}$. If $p_{\ell}>1-\alpha$, then $\epsilon_{\ell+1}=(\epsilon_{\ell}+\epsilon_{LB})/2$ and $\epsilon_{UB}=\epsilon_{\ell}$; if $p_{\ell}<1-\alpha$, then $\epsilon_{\ell+1}=(\epsilon_{UB}+\epsilon_{\ell})/2$, or $\epsilon_{\ell+1}=2\epsilon_{0}$ if $\epsilon_{UB}=\infty$, and $\epsilon_{LB}=\epsilon_{\ell}$. We now solve the problem again for $\epsilon_{\ell+1}$ and compute $p_{\ell+1}$. We continue in this fashion until the estimated probability, $p_{\ell}$, is within a threshold of the target probability $1-\alpha$.

One way to pick the initial $\epsilon_{0}$ for the binary search is as follows. First, compute $x_{0}^{*}$ as the optimal solution of the “robust” formulation of Eq. 1,

[TABLE]

Then, evaluate the standard deviation of the sample $\{C(x_{0}^{*},\xi_{1}),\ldots,C(x_{0}^{*},\xi_{N})\}$, and set $\epsilon_{0}$ to be twice the standard deviation. This way of selecting $\epsilon_{0}$ usually renders conservative solutions; however, because we have observed that the optimization algorithm converges more quickly for large values of $\epsilon$, we prefer to start with a large value of $\epsilon$ that is decreased afterwards. This choice of $\epsilon_{0}$ is motivated by the attempt to minimize the mean square error in Eq. 9, for a fixed $x$. In **[1]**, the authors show that the best choice of $\epsilon$ is proportional to its standard deviation if everything is assumed to be normal. Thus, we define $\epsilon_{0}$ as a multiple of the standard deviation evaluated at $x_{0}^{*}$.

To estimate the true probability of constraint violation, we obtain an out-of-sample approximation using a sample with a large number of scenarios $N^{\prime}$ and compute the quantity $\sum_{i=1}^{N^{\prime}}\mathbbm{1}(C(x,\xi_{i})\leq 0)/N^{\prime}$ as an approximation of $\operatorname{\mathbb{P}}(C(x,\xi)\leq 0)$.

To speed up the binary search, we use the solution $x^{*}_{\epsilon_{k}}$, with its multipliers $\lambda^{*}_{\epsilon_{k}}$ and $\bar{\mu}^{*}_{\epsilon_{k}}$, to initialize the algorithm for $\epsilon_{k+1}$, i.e., $x_{0}=x^{*}_{\epsilon_{k}}$, $\lambda_{0}=\lambda^{*}_{\epsilon_{k}}$, and $\bar{\mu}_{0}=\bar{\mu}^{*}_{\epsilon_{k}}$ to solve problem Eq. 20 for $\epsilon_{k+1}$.

For the experiments in the following section, we terminate the binary search if $|\operatorname{\mathbb{P}}(C(x_{\ell}^{*},\xi)-(1-\alpha)|\leq 10^{-4}$ or the number of bisections exceeds 10.

6.2 Single Chance Constraints

In this section, we present the outcome of two experiments using a nonlinear solver to solve our reformulation for a SCCP. Both examples are solved by Knitro 10.1.2. To compute $Q_{\epsilon}$, we use the Matlab function fzero to find the root of the nonlinear function Eq. 7 from the initial point $C_{[M]}(x)$. We also provide Knitro both the gradient and Hessian of $Q_{\epsilon}$ as given in Eq. 17 and Eq. 18.

6.2.1 A non-convex example

We first investigate the performance of the NLP reformulation for the non-convex problem

[TABLE]

where $c(x,\xi)=0.25x^{4}-1/3x^{3}-x^{2}+0.2x-19.5+\xi_{1}x+\xi_{2},$ and $\xi_{1}\sim N(0,3)$ and $\xi_{2}\sim N(0,144)$ are independent random variables. This problem is similar to the one presented in **[12]**.

In this experiment we use a sample of size $N=1,000$ and two different smoothing parameters: $\epsilon=1$ and $\epsilon=0.1$. We used 10 different feasible initial points for each $\epsilon$, given by $x_{0}=(t,2.5)$, for $t$ equally spaced in the interval $[-1.5,2.5]$. The solving time for all starting points and both $\epsilon$ values is less than 0.0867 seconds. In Fig. 3 we observe the results of all the different runs. In this figure, the solid line represents the boundary of the empirical feasible region, while the dotted line represents the boundary of the smooth approximation to the feasible region. All the points above the solid and dotted line are feasible for the empirical and smooth approximations, respectively. We see that the smaller value of $\epsilon$ translates into a smooth curve that follows the empirical quantile more closely, while the larger value renders a “less noisy” function.

Figure 3* shows that our approach is designed to find local minima instead of global minima, and that the quality of the solution not only depends on the starting point but on the choice of $\epsilon$. This can be seen from the points in Fig. 3. Notice that the true problem seems to have only two local minimizers, one around $x=1.8$ and another around $x=-0.8$. However, when $\epsilon=0.1$, the solver gets attracted to spurious local minima for three different starting points.*

It is also interesting to notice that the algorithm takes longer to find a solution the smaller $\epsilon$ is. For $\epsilon=1$, the number of iterations ranged between 6 and 10 with an average of 8.2; however, for $\epsilon=0.1$, the number of iterations ranges between 7 and 14 with an average of 11.6. This exemplifies some of the claims made in section 2.1 and emphasizes the importance of choosing $\epsilon$ appropriately.

6.2.2 Reinsurance portfolio example

We consider a Value-at-Risk (VaR) minimization problem for reinsurance portfolios. The decision variable $x_{j}$, for $j=1,\ldots,n$, represents the fraction of the risk contract $j$ that the reinsurer is willing to take. We would like to choose the optimal fraction of the contracts such that the portfolio premium satisfies a pre-established minimum value and that the $0.95$-VaR is minimized. Each contract is subject to a stochastic loss. This data is simulated to represent risk and returns of businesses. The characteristics of these loss matrices include non-negativity, sparsity, positive skewness, and high kurtosis (this problem was constructed as in **[14]**).

Notice that the $0.95$-VaR is defined the same way as the $0.95$-quantile, so we can therefore use $Q^{0.95}$ as the objective function. Thus, the problem can be modeled as

[TABLE]

where the random variable $L_{j}$ represents the losses of contract $j$, $c_{j}$ are the market premiums of the risk contracts, and $c^{*}=0.1\sum_{j=1}^{n}c_{j}$ is the minimum level of portfolio premium. The chance-constrained reformulation of the problem is

[TABLE]

For this example, we considered $25$ contracts ($n=25$) for different sample sizes $N$, and we ran 10 replications for each $N$.

In Table 2, we present the statistics of the Knitro runs using formulation Eq. 8, the quantile-constrained NLP approximation of the problem. The times reported include all NLP solves executed during the bisection search to tune the smoothing parameter $\epsilon$ described in section 6.1. We see that the solution time grows modestly with respect to the sample size, and that the number of iterations are roughly the same even when the sample size increases. The last column displays the average value of the tuned smoothing parameter. Consistent with our theoretical analysis in Section 3 and the results in **[2]**, we observe that the value of $\epsilon$ that yields a solution that is “just feasible” decreases as the sample size increases.

We compare our approach to the empirical probability formulation $\widetilde{F}^{N}$. This formulation is solved via a mixed-integer linear programming (MIP) formulation **[38]**, using Gurobi as the MIP solver with a time limit of 10 minutes. Mixing inequalities **[28]** were not used, because the required preprocessing took already more than 10 min for $N=2,000$.

*For the solutions obtained by either method, we computed an out-of-sample approximation of the true quantile, $Q^{0.95}(L^{T}x)$, via the empirical $0.95$-quantile using a sample of *$N^{\prime}=$100,000 observations. The results of the experiments are shown in Fig. 4. We see that, for $N=200$ and $N=5,000$, the variance of the solutions obtained by the NLP solver is considerably smaller than the variance of the solutions obtained using the MIP formulation. The spread observed on the quality of the solutions from the MIP for $N=5,000$ may be due to the 10 minute limit, since the incumbent returned by the solver is not proved to be optimal for the sample approximation for this sample size. We also observe that the solutions obtained by the NLP formulation are consistently better in all cases and are improving as the sample size grows.

We now compare the performance of the NLP solver when constraint Eq. 16b is substituted by $F_{\epsilon}(0;x)=\frac{1}{N}\sum_{i=1}^{N}\Gamma_{\epsilon}(c(x,\xi_{i}))\geq 1-\alpha$. We use this example to evaluate the advantage of the quantile formulation of the constraint compared to the probability formulation. We did these experiments under the same settings as above. The solutions that the solver found with the probability formulation are the same as for the quantile formulation (within numerical accuracy), except for two cases for $N=2,000$ and one for $N=5,000$ in which Knitro’s iterates diverge. The number of iterations and the time taken by the solver is considerably larger compared to the quantile formulation as can be seen in Fig. 5. Furthermore, while the quantile formulation requires roughly the same number of iterations for different number of scenarios, the iteration count for the probabilistic constraint grows with the sample size. Therefore, although finding the root of Eq. 7 is more expensive than using $F_{\epsilon}(0;x)$ directly, working with the quantile is preferable in terms of the performance of the solver.

6.3 Joint Chance Constraints

The last example we consider is a continuous probabilistic multi-dimensional knapsack problem

[TABLE]

The problem consists of a set of $n$ items, each with an associated profit $p_{j}$, and different “weights” for each item (for example, both volume and weight, where the volume and weight of each item are not necessarily related). Here, $W_{ij}$ represents the random $i$th “weight” of item $j$, and $w_{i}$ is the maximum “weight” limit for constraint $i$ ($w_{i}$ is a deterministic quantity). The decision variable $x_{j}\in[0,1]$ represents the fraction of an item that is to be included. The goal is to select the optimal fraction each item in order to maximize the profit, while satisfying all the “weight” constraints simultaneously at least $95\%$ of the time. The number of items considered for this experiment was $n=20$ with $m=10$ constraints. The randomness comes from the weights $W_{ij}$, whose distribution is described next.

We implemented Algorithm 1 in Matlab, using CPLEX 12.6.3 to solve the QP Eq. 27. We selected the following parameters for the algorithm: $\pi=10$, $\hat{\Delta}=10^{6}$, $\Delta_{0}=1$, $\eta=10^{-8}$, $\tau_{1}=1/2$ and $\tau_{2}=2$. We defined $H^{k}$ as Eq. 36, and terminate Algorithm 1 when $\|\nabla\overline{\mathcal{L}}(x_{k},\nu,\lambda)\|_{\infty}\leq 10^{-6}$, $g(x_{k})\leq 10^{-6}\boldsymbol{1}_{m}$, and $Q_{\epsilon}(C^{N}(x_{k}))<10^{-6}$.

We generate data for Eq. 42 according to the knapsack instance mk-20-10 from **[41]**, which is based on an available deterministic knapsack instance in **[3]**. Each item $j$ is available or not according to a Bernoulli distribution. An item that is not available has weight zero in all rows. As a result, the item weights in different rows are not independent. If an item is available, then its weight in each row is normally distributed with mean equal to its weight in the deterministic instance and standard deviation equal 0.1 times the mean. For each sample size $N$, we ran 10 replications.

The results of our proposed algorithm are presented in Table 3 under the column heading ‘NLP’. The computation times include the total time for the binary search algorithm to find a tuned value of $\epsilon$. For the solutions obtained, we computed an out-of-sample approximation of $\operatorname{\mathbb{P}}(C(x,\xi))$ using the empirical cdf with $N^{\prime}=10^{6}$ scenarios. First, we can see that the time it takes the algorithm to obtain a solution grows modestly with the sample size $N$. Second, the best solution obtained from solving the problem with 10 different samples of size 100 has an objective value nearly identical to the best solution obtained using samples of size 1,000. This indicates that solving the approximation with multiple small samples may be more effective at obtaining good solutions than solving the approximation with a single large sample. Finally, the last row shows that the average value of the tuned smoothing parameter decreases with increased sample size, similar to our observation in Table 2

We also compared our solutions to the solutions obtained by solving the empirical probability formulation $\widetilde{F}^{N}$ via an MIP formulation. In order to obtain comparable solutions between the MIP and the NLP formulations, we tuned the risk level $\alpha$ of the MIP formulation to ensure that the solutions obtained were all feasible. This is similar to tuning the $\epsilon$ parameter in the NLP formulation. We limited the overall tuning time for the $\alpha$ parameter in the MIP to ten minutes. For each $\alpha$ trial of the MIP, we gave Gurobi 0.5, 0.5, 1, 5, and 10 minutes for sample size 100, 500, 1,000, 5,000, and 10,000, respectively. These time limits were chosen to give the MIP algorithm significantly more time than the NLP approach. The results are given in Table 3 under the column heading ‘MIP’. It can be seen that the MIP solutions have generally lower optimal values and larger variability, the same behavior observed in section 6.2.2. On the other hand, solutions obtained from Algorithm 1 consistently present larger objective function values, less variability, and are solved well under 3 minutes for most sample sizes (except $N=10,000$). It is possible that giving the MIP algorithm more time could lead to better solutions, but these results clearly indicate that our NLP approach is more effective at finding good solutions quickly. Furthermore, the results show that the binary search for $\epsilon$ can attain the desired risk-level more accurately than changing the $\alpha$ parameter, given the time limit.

7 Summary

In this paper we developed a technique to solve nonlinear problems with probabilistic constraints. We first presented an approximation for single chance constraints with the following characteristics: (1) it can be represented as a single smooth constraint, (2) it does not require information about the specific distribution of the random variables or information about the moments since it solely based on sampling; it only needs to satisfy continuity of the distribution, (3) it does not restrict the constraints to be convex, (4) it can be handled by standard continuous nonlinear programming techniques and solvers, and (5) it avoids the combinatorial complexity of a branch and bound search, and scales well with the number of variables and samples.

We also presented an extension of the SCC formulation that is able to handle joint chance constraints. We constructed a specialized trust-region solver for JCCPs utilizing a Hessian matrix that shows fast local convergence in our experiments.