I want to solve an optimization model that contains a constraint like $$ \Pr[F(x,\xi)\leq0]\geq1-\varepsilon $$ where $x$ are my decision variables, $\xi$ is a random vector, and $\varepsilon\in(0,1)$ is some tolerance that measures how certain I want to be to satisfy my constraint. Is there a general way to solve such a problem?
(Disclaimer: I'm asking this because the Nemirovski and Shapiro paper was hard for me to read, and I wanted a brief, accessible summary).
In general, chance constraints are difficult to handle computationally. Given a fixed $x$, is it often very difficult to even compute the value of $\Pr[F(x,\xi)\leq0]$, making the incorporation of such a constraint into an optimization model quite troublesome.
One possible approach is outlined in the 2006 paper by Nemirovksi and Shapiro entitled Convex Approximations of Chance Constrained Programs. In this paper, the authors propose techniques for creating conservative approximations of chance constraints (i.e., approximations with smaller feasible regions). See this paper by Shabbir Ahmed for convex relaxations.
In what follows, we'll assume that for all $\xi$, the function $F(\cdot,\xi)$ is convex in its first argument (if this is not the case, you'll have to do some hard thinking). We'll also assume $F$ is real-valued (not vector-valued---i.e. not a joint chance constraint). For the remainder, we'll deal with the complement of the given constraint, that is, we'll deal with $\Pr[F(x,\xi)>0]\leq\varepsilon$. It is well known that the probability of an event is equal to the expectation of an indicator function on that event. That is, $$ \Pr[F(x,\xi)>0]=\mathbb{E}[\mathbf{1}_{\{F(x,\xi)>0\}}]. $$ The BIG IDEA behind this approximation scheme is to replace the discontinuous, non-convex indicator function $\mathbf{1}_{\{\cdot\}}$ with some nicer function---let's call it $\psi(\cdot)$. Ideally, we want $\psi$ to be continuous, convex (for easier computing), and be greater than (or equal to) the indicator function (so that our approximation is conservative).
Now, because $\psi(t)\geq\mathbf{1}_{\{t>0\}}$, we have that $$ \Pr[F(x,\xi)>0]=\mathbb{E}[\mathbf{1}_{\{F(x,\xi)>0\}}]\leq\mathbb{E}[\psi(F(x,\xi))]. $$ Hence, if we enforce $\mathbb{E}[\psi(F(x,\xi))]\leq\varepsilon$, this implies $\Pr[F(x,\xi)>0]\leq\varepsilon$ (just like we want).
So how do we pick $\psi$, the approximation to the indicator function? Here are three popular choices:
In fact, you can do a little bit better than this approximation. Nemirovski and Shapiro show that you can actually enforce $$ t\mathbb{E}[\psi(\tfrac{1}{t}F(x,\xi))]\leq{t}\varepsilon\quad\leftarrow\text{If you only remember one thing...} $$ for any $t>0$, and the approximation is still valid. This function on the left-hand side is actually convex in $t$ (as long as $\psi$ is convex)---this is called a perspective function (see, e.g., Section 3.2.6 of Boyd and Vandenberghe's Convex Optimization). So we can just throw in a new variable $t\geq0$ and enforce the above constraint, and get a tighter approximation.
Toy Example
Suppose we have the constraint $$ \Pr[\xi_1x_1+\xi_2x_2\leq2]\geq0.9. $$ where $x_1$ and $x_2$ are non-negative decision variables, $\xi_1\sim\mathrm{Poisson}(1)$, and $\xi_2\sim\mathrm{Poisson}(1/2)$. The true chance constraint can be evaluated using the law of total probability (see code below). For the conservative approximation, we take $F(x,\xi)=\xi_1x_1+\xi_2x_2-2$, $\varepsilon=0.1$ and $\psi(t)=\exp(t)$. As noted above, this choice of $\psi$ is equivalent to the MGF of the random variable $\xi_1x_1+\xi_2x_2-2$, which is given by $$ \exp\big[\lambda_1(e^{x_1/t}-1)+\lambda_2(e^{x_2/t}-1)-2/t\big] $$ where $\lambda_1=1$ and $\lambda_2=1/2$. Hence, for any $t>0$, the following is a conservative approximation: $$ t\cdot\exp\big[(e^{x_1/t}-1)+\tfrac{1}{2}(e^{x_2/t}-1)-2/t\big]\leq{t\varepsilon} $$ An example of the true feasible region (light gray), and the approximation for $t=0.4$ (dark gray) are shown here:
Python code to produce this plot:
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def main():
l1, l2 = 1.0, 0.5
tau = 2. # Chance constraint rhs
eps = 0.1 # Chance risk tolerance
tee = 0.4 # Conservatism a parameter
num_points = 500
X = np.linspace(0, 1, num_points)
X, Y = np.meshgrid(X, X)
approx = tee * np.exp(
l1 * (np.exp(X / tee) - 1) + l2 * (np.exp(Y / tee) - 1) - (tau / tee)
)
exact = np.vectorize(evaluate_poisson_cc)(X, Y, l1, l2, tau)
# Plot the results
ax = plt.gca()
ax.contourf(X, Y, exact, levels=[1. - eps, 1.], colors="k", alpha=0.2)
ax.contourf(X, Y, approx, levels=[0., tee * eps], colors="k", alpha=0.4)
ax.set_xlabel("$x_1$"); ax.set_ylabel("$x_2$")
ax.set_aspect("equal")
plt.show()
def evaluate_poisson_cc(x1, x2, l1, l2, tau, tol=1e-10):
""" Evaluate P[x1 * Z1 + x2 * Z2 <= tau] where Zi ~ Poisson(li)
"""
assert tau > tol
if np.abs(x1) < tol and np.abs(x2) < tol:
return 1.
if np.abs(x1) < tol:
return stats.poisson.cdf(tau / x2, l2)
if np.abs(x2) < tol:
return stats.poisson.cdf(tau / x1, l1)
max_value = np.floor(tau / x2)
n_values = np.arange(max_value + 1)
v1 = stats.poisson.pmf(n_values, l2)
v2 = stats.poisson.cdf((tau - n_values * x2) / x1, l1)
return np.dot(v1, v2)
if __name__ == "__main__":
main()
In addition to David M's excellent answer, another approach which sometimes can be useful is to use integer programming. Specifically, if your distribution of $\xi$ is discrete (e.g. taking a sample-average-approximation of uncertainty is a reasonable approach) and you are willing to make a big-M assumption on $F(x, \xi)$ then another way to model a chance constraint is using integer programming.
In particular, if we have $L$ realizations of $\xi_i$ total, and they all occur with equal probability, then we can let $z_i=1$ if $F(x, \xi_i) \leq 0$ and $0$ otherwise, and model the chance constraint exactly using the following big-M constraints: $$ F(x, \xi_i) \leq M(1-z_i), \forall i \in \{1, .., L\},\\ \frac{1}{L}\sum_i z_i \geq 1-\epsilon. $$ This approach cannot be applied to a general continuous distribution, because it relies on the support of $\xi$ being discrete. However, there are some fairly general conditions under which the optimal solutions $x^*_L$ from sample-average-approximations of size $L$ converge almost-surely to the optimal solution $x^*$ for the full continuous distribution of $\xi$ as $L \mapsto \infty$, and therefore discretizing a continuous distribution $\xi$ and applying a big-M approach to the discretized distribution is often a reasonable idea.
I believe that some other approaches have since been shown to be faster (see here for instance), but the integer programming approach seems worth knowing.
Further reading: chapter 5.7 of this book, and this paper.
