Representation of chance-constraints with strong asymptotic guarantees

TL;DR

This paper develops a method to approximate chance-constrained feasible sets using polynomial outer and inner approximations with strong asymptotic guarantees, ensuring convergence in measure as polynomial degree increases.

Contribution

It introduces a sequence of polynomial-based outer and inner approximations for chance-constraints with proven convergence guarantees in measure.

Findings

Outer approximations converge to the true feasible set in measure.

Inner approximations also converge with similar guarantees.

The approach uses semidefinite programming to compute the approximations.

Abstract

Given $\epsilon \in (0,1)$, a probability measure $\mu$ on $\Omega\subset\mathbb{R}^p$ and a semi-algebraic set $K\subset X\times\Omega$, we consider the feasible set $X^*_\epsilon=\{x\in X:{\rm Prob}[(x,\omega)\in K]\geq 1-\epsilon\}$ associated with a chance-constraint. We provide a sequence of outer approximations $X^d_\epsilon=\{x\in X: h_d(x)\geq0\}$, $d\in\mathbb{N}$, where $h_d$ is a polynomial of degree $d$ whose vector of coefficients is an optimal solution of a semidefinite program. The size of the latter increases with the degree $d$. We also obtain the strong and highly desirable asymptotic guarantee that $\lambda(X^d_\epsilon\setminus X^*_\epsilon)\to0$ as $d$ increases, where $\lambda$ is the Lebesgue measure on $X$. Inner approximations with same guarantees are also obtained.

Figures

Tables

Table 1. TABLE I: The effect of Stokes constraints. n = 2 𝑛 2 n=2

	${\rho }_3$	${\rho }_3^{\prime }$	${\rho }_4$	${\rho }_4^{\prime }$
a =1.4	5.71	3.55	5.38	3.27
a= 1.3	5.38	3.41	5.04	3.21
a= 1.2	5.02	3.31	4.70	3.18
a= 1.1	4.56	3.25	4.32	3.16
a= 1.0	3.91	3.20	3.87	3.15

Table 2. TABLE II : Example III.5 ; n=1: with and without Stokes

	d=4	d=4 (Stokes)	d=8	d=8 (Stokes)
$ϵ=0.75$	13%	3.8%	6.3%	0.7%
$ϵ=0.50$	16.6%	2.1%	9.8%	0.2%
$ϵ=0.25$	26.6%	6.5%	18.5%	1.0%
$ϵ=0.00$	44.7%	22.7%	31.2%	8.2%

Table 3. TABLE III : Example III.6 , n=1; with and without Stokes

	d=7	d=7 (Stokes)	d=10	d=10 (Stokes)
$ϵ=0.7$	2.3%	0.4%	2.3%	0.2%
$ϵ=0.4$	30%	10%	26%	5.5%
$ϵ=0.1$	52%	24%	46%	16%

Equations

$$\mathbf{K}=\{(\mathbf{x},\omega ):g_j(\mathbf{x},\omega )\ge 0,j=1,\dots ,m\},$$

$${\mathbf{K}}_{\mathbf{x}}:=\{\omega \in \mathbf{\Omega }:(\mathbf{x},\omega )\in \mathbf{K}\},\mathbf{x}\in \mathbf{X}.$$

$${\mathbf{X}}_{ϵ}^{\ast }$$

$${\mathbf{X}}_{ϵ}^d:=\{\mathbf{x}\in \mathbf{X}:h_d(\mathbf{x})\ge 0\},d\in \mathbb{N},$$

$$\min \{f(\mathbf{x}):\mathbf{x}\in C;\mathrm{Prob}((\mathbf{x},\omega )\in \mathbf{K})\ge 1-ϵ\}.$$

$$\min \{f(\mathbf{x}):\mathbf{x}\in C;h_d(\mathbf{x})\ge 0\},$$

$${\mathbf{x}}^{\ast }=\arg \max \{\mu ({\mathbf{K}}_{\mathbf{x}}):\mathbf{x}\in \mathbf{X}\}.$$

$$\underset{d\to \infty }{\lim }\lambda ({\mathbf{X}}_{ϵ}^d\setminus {\mathbf{X}}_{ϵ}^{\ast })=0,$$

$$\mathbf{x}\mapsto f(\mathbf{x})=\sum _{\alpha \in {\mathbb{N}}^n}{f}_{\alpha }{\mathbf{x}}^{\alpha }.$$

$$f(=\sum _{\alpha }{f}_{\alpha }{\mathbf{x}}^{\alpha })\mapsto L_{\mathbf{y}}(f):=\sum _{\alpha }{f}_{\alpha }{y}_{\alpha }.$$

$$M_d(\mathbf{y})(\alpha ,\beta ):=L_{\mathbf{y}}({\mathbf{x}}^{\alpha +\beta })=y_{\alpha +\beta },\alpha ,\beta \in {\mathbb{N}}_d^n.$$

$$M_d(g\mathbf{y})(\alpha ,\beta )$$

$$\lambda (\mathbf{K})=\underset{\varphi \in M(\mathbf{K})}{\sup }\{\varphi (\mathbf{K}):\varphi \le \lambda \},$$

$$\rho$$

$$=$$

$$\begin{array}[]{rl}\rho_{d}\,=\,\displaystyle\sup_{\mathbf{y},\mathbf{z}}&\{\,L_{\mathbf{y}}(1):y_{\alpha}+z_{\alpha}=\lambda_{\alpha},\quad\forall\alpha\in\mathbb{N}^{n}_{2d}\\ &\mathbf{M}_{d}(\mathbf{y}),\mathbf{M}_{d}(\mathbf{z})\succeq 0\\ &\mathbf{M}_{d-d_{j}}(g_{j}\,\mathbf{y})\succeq 0,\quad j=1,\ldots,m\}.\end{array}$$

$${\int }_{\mathbf{K}}\underset{\mathbf{x}\mapsto {\theta }_{\alpha }^i(\mathbf{x})}{\underbrace{\frac{\partial }{\partial x_i}({\mathbf{x}}^{\alpha }h(\mathbf{x}))}}d\lambda (\mathbf{x})=0,i=1,\dots ,n,$$

$${\int }_{\mathbf{K}}{\theta }_{\alpha }^i(\mathbf{x})d{\varphi }^{\ast }(\mathbf{x})=0,\forall \alpha \in {\mathbb{N}}^n;i=1,\dots ,n.$$

$$L_{\mathbf{y}}({\theta }_{\alpha }^i)=0,\forall \alpha \in {\mathbb{N}}_{2d-\mathrm{deg}(h)}^n;i=1,\dots ,n.$$

$$\theta_{\alpha}(\mathbf{x})\,=\,\left\{\begin{array}[]{ll}(n+d_{f})\mathbf{x}^{\alpha}\,h+\langle\mathbf{x},\nabla(\mathbf{x}^{\alpha}\,h)\rangle&\mbox{for $fd\lambda$}\\ \mathbf{x}^{\alpha}\,h(n+d_{f}f)+\langle\mathbf{x},\nabla(\mathbf{x}^{\alpha}\,h)\rangle&\mbox{for $\exp(f)d\lambda$}.\end{array}\right.$$

$$\lambda \otimes \mu (A\times B)=\lambda (A)\mu (B),\forall A\in \mathcal{B}(\mathbf{X}),B\in \mathcal{B}(\mathbf{\Omega }).$$

$$\rho =\underset{\varphi \in M(\mathbf{K})}{\sup }\{\varphi (\mathbf{K}):\varphi \le \lambda \otimes \mu \}.$$

$$d{\varphi }^{\ast }((\mathbf{x},\omega ))=1_{\mathbf{K}}((\mathbf{x},\omega ))d\lambda \otimes \mu ((\mathbf{x},\omega )),$$

$$\rho$$

$${\int }_{\mathbf{X}\times \mathbf{\Omega }}{1}_{\mathbf{K}}((\mathbf{x},\omega ))\lambda \otimes \mu (d(\mathbf{x},\omega ))=$$

$${\int }_{\mathbf{X}}\underset{\mu ({\mathbf{K}}_{\mathbf{x}})}{\underbrace{\left({\int }_{\mathbf{\Omega }}{1}_{\mathbf{K}}((\mathbf{x},\omega ))\mu (d\omega )\right)}}\lambda (d\mathbf{x})={\int }_{\mathbf{X}}\mu ({\mathbf{K}}_{\mathbf{x}})\lambda (d\mathbf{x})$$

$$\begin{array}[]{rl}\rho_{d}\,=\,\displaystyle\sup_{\mathbf{y},\mathbf{z}}&\{\,y_{0}:\\ \mbox{s.t.}&y_{\alpha,\beta}+z_{\alpha,\beta}\,=\,\lambda_{\alpha}\cdot\mu_{\beta},\quad(\alpha,\beta)\in\mathbb{N}^{n+p}_{2d}\\ &\mathbf{M}_{d}(\mathbf{y}),\,\mathbf{M}_{d}(\mathbf{z})\succeq 0\\ &\mathbf{M}_{d-d_{j}}(g_{j}\,\mathbf{y})\,\succeq\,0,\quad j=1,\ldots,m\},\end{array}$$

$$\begin{array}[]{rl}\rho^{*}_{d}\,=\,\displaystyle\inf_{p\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]_{2d}}&\{\displaystyle\int_{\mathbf{X}\times\mathbf{\Omega}}p(\mathbf{x},\mathbf{\omega})\,\lambda\otimes\mu(d(\mathbf{x},\mathbf{\omega}))\\ \mbox{s.t.}&p(\mathbf{x},\mathbf{\omega})\,\geq\,1,\quad\forall(\mathbf{x},\mathbf{\omega})\in\mathbf{K}\\ &\mbox{$p$ is SOS}\,\}.\end{array}$$

$${\rho }_d={\rho }_d^{\ast }={\int }_{\mathbf{X}\times \mathbf{\Omega }}{p}_d^{\ast }(\mathbf{x},\omega )\lambda \otimes \mu (d(\mathbf{x},\omega )).$$

$$\mathbf{x}\mapsto h_d^{\ast }(\mathbf{x}):={\int }_{\mathbf{\Omega }}{p}_d^{\ast }(\mathbf{x},\omega )\mu (d\omega ),\mathbf{x}\in {\mathbb{R}}^n.$$

$${\rho }_d={\int }_{\mathbf{X}}{h}_d^{\ast }(\mathbf{x})\lambda (d\mathbf{x})\to \rho ={\int }_{\mathbf{X}}\mu ({\mathbf{K}}_{\mathbf{x}})\lambda (d\mathbf{x})$$

Full text

I Introduction

We consider the following general framework for decision under uncertainty : Let $\mathbf{x}\in\mathbf{X}\subset\mathbb{R}^{n}$ be a decision variable while $\mathbf{\omega}\in\mathbb{R}^{p}$ is a disturbance (or noise) parameter whose distribution $\mu$ (with support $\mathbf{\Omega}\subset\mathbb{R}^{p}$) is known, i.e., its list of moments $\mu_{\beta}:=\int_{\mathbf{\Omega}}\mathbf{\omega}^{\beta}\,d\mu(\mathbf{\omega})$, $\beta\in\mathbb{N}^{p}$, is available in closed form or numerically.

Both $\mathbf{x}$ and $\mathbf{\omega}$ are linked by constraints of the form $(\mathbf{x},\mathbf{\omega})\in\mathbf{K}\subset\mathbf{X}\times\mathbf{\Omega}$, where

[TABLE]

for some polynomials $(g_{j})\subset\mathbb{R}[\mathbf{x},\mathbf{\omega}]$, that is, $\mathbf{K}$ is a basic semi-algebraic set.

Next, for each fixed $\mathbf{x}\in\mathbf{X}$, let $\mathbf{K}_{\mathbf{x}}\subset\mathbf{\Omega}$ be the (possibly empty) set defined by:

[TABLE]

Let $\epsilon\in(0,1)$ be fixed. The goal of this paper is to provide tight approximations of the set

[TABLE]

in the form:

[TABLE]

where $h_{d}$ is a polynomial of degree at most $d$.

Such approximations are particularly useful for optimization and control problems with chance-constraints; for instance problems of the form:

[TABLE]

Indeed one then replaces problem (5) with

[TABLE]

where the uncertain parameter $\mathbf{\omega}$ has disappeared. So if $C$ is a basic semi-algebraic set then (6) is a standard polynomial optimization problem. Of course the resulting decision problem (6) may still be hard to solve because the sets $\mathbf{X}^{*}_{\epsilon}$ and $\mathbf{X}^{d}_{\epsilon}$ are not convex in general. But this may be the price to pay for avoiding a too conservative formulation of the problem. The interested reader is referred to Henrion [3], Prékopa [14] and Shapiro [15] for a general overview of chance constraints in optimization and to Calafiore and Dabbene [2], Jasour et al. [7] and Li et al. [13] in control (and the references therein).

However, in the formulation (6) one has got rid of the disturbance parameter $\mathbf{\omega}$, and so one may apply the arsenal of Non Linear Programming algorithms to get a local minimizer of (6). If $n$ is not too large or if some sparsity is present in problem (6) one may even run a hierarchy of semidefinite relaxations to approximate its global optimal value. For the latter approach the interested reader is referred to [10] and for a discussion on this approach to various control problems with chance constraints we refer to the recent paper of Jasour et al. [7] and the references therein.

In Jasour et al. [7] the authors have considered some control problems with chance constraints. They have provided an elegant formulation and a numerical scheme for solving the related problem of computing

[TABLE]

This problem is posed as an infinite-dimensional LP problem in an appropriate space of measures, that is, a Generalized Moment Problem (GMP) as described in Lasserre [10]. Then to obtain $\mathbf{x}^{*}$ they solve a hierarchy of semidefinite relaxations, which is the moment-SOS approach for solving the GMP. This GMP formulation has the particular and typical feature of including a constraint of domination $\phi\leq\psi$ between two measures $\phi$ and $\psi$. Such domination constraints are particularly powerful and have been already used in a variety of different contexts. See for instance Henrion et al. [5] for approximating the Lebesgue volume of a compact semi-algebraic set, Lasserre [12] for computing Gaussian measures of semi-algebraic set, Lasserre [11] for “approximating” the Lebesgue decomposition of a measure with respect to another one. It has been used by Henrion and Korda [6] for approximating regions of attraction, by Korda et al. [9] for approximating maximum controlled invariant sets, and more recently in Jasour and Lagoa [8] for a unifying treatment of some control problems.

Contribution

The approach that we propose for determining the set $\mathbf{X}^{d}_{\epsilon}$ defined in (4) is very similar in spirit to that in [5] and [7] and can be viewed as an additional illustration of the versatility of the GMP and the moment-SOS approach in control related problems. Indeed we also define an infinite-dimensional LP problem $\mathbf{P}$ in an appropriate space of measures and a sequence of semidefinite relaxations $(\mathbf{P}_{d})_{d\in\mathbb{N}}$ of $\mathbf{P}$, whose associated monotone sequence of optimal values $(\rho_{d})_{d\in\mathbb{N}}$ converges to the optimal value $\rho$ of $\mathbf{P}$. An optimal solution of the dual of $(\mathbf{P}_{d})$ allows to obtain a polynomial $h_{d}$ of degree $2d$ whose super-level set $\{\mathbf{x}:h_{d}(\mathbf{x})\geq 0\}$ is precisely the desired approximation $\mathbf{X}^{d}_{\epsilon}$ of $\mathbf{X}^{*}_{\epsilon}$ in (4); in fact the sets $(\mathbf{X}^{d}_{\epsilon})_{d\in\mathbb{N}}$ provide a sequence of outer approximations of $\mathbf{X}^{*}_{\epsilon}$. We also provide the strong asymptotic guarantee that

[TABLE]

where $\lambda$ is the Lebesgue measure on $\mathbf{X}$, which to the best of our knowledge is the first result of this kind at this level of generality. (The same methodology applied to chance-constraints of the form ${\rm Prob}((\mathbf{x},\mathbf{\omega})\in\mathbf{K})<\epsilon$ would provide a sequence of inner approximations of the set $\{\mathbf{x}\in\mathbf{X}:{\rm Prob}((\mathbf{x},\mathbf{\omega})\in\mathbf{K})<\epsilon\}$.)

Another contribution is to include a technique to accelerate the convergence $\rho_{d}\to\rho$ which otherwise can be too slow. This technique is different from the one used in [5] for the related problem of computing the volume of a semi-algebraic set, and has the nice feature or preserving the monotonicity of the convergence of $\rho_{d}\to\rho$. It can be applied whenever $d\mu$ is the Lebesgue measure $d\mathbf{\omega}$ on $\mathbf{\Omega}$ or $d\mu=q(\mathbf{\omega})\,d\mathbf{\omega}$ or $d\mu=\exp(-q(\mathbf{\omega}))d\mathbf{\omega}$ for some homogeneous nonnegative polynomial $q$.

At last but not least, in principle we can also treat the case where the support $\mathbf{\Omega}$ of $\mu$ and the set $\mathbf{K}$ are not compact, which includes the important case where $\mu$ is the normal distribution. We briefly explain what are the (technical) arguments which allow to extend the method to the non compact case.

Of course this methodology is computationally expensive and so far limited to relatively small size problems (but after all the problem is very hard). An interesting issue not discussed here is to investigate whether sparsity patterns can be exploited to handle problems with larger size.

II Notation, definitions and preliminary results

II-A Notation and definitions

Let $\mathbb{R}[\mathbf{x}]$ be the ring of polynomials in the variables $\mathbf{x}=(x_{1},\ldots,x_{n})$ and let $\mathbb{R}[\mathbf{x}]_{d}$ be the vector space of polynomials of degree at most $d$ (whose dimension is $s(d):={n+d\choose n}$). For every $d\in\mathbb{N}$, let $\mathbb{N}^{n}_{d}:=\{\alpha\in\mathbb{N}^{n}:|\alpha|\,(=\sum_{i}\alpha_{i})\leq d\}$, and let $\mathbf{v}_{d}(\mathbf{x})=(\mathbf{x}^{\alpha})$, $\alpha\in\mathbb{N}^{n}_{d}$, be the vector of monomials of the canonical basis $(\mathbf{x}^{\alpha})$ of $\mathbb{R}[\mathbf{x}]_{d}$. A polynomial $f\in\mathbb{R}[\mathbf{x}]_{d}$ is written

[TABLE]

Given a closed set $\mathcal{X}\subset\mathbb{R}^{n}$, denote by $M(\mathcal{X})$ the space of finite Borel measures on $\mathcal{X}$ and by $\mathcal{P}(\mathcal{X})$ the convex cone of polynomials that are nonnegative on $\mathcal{X}$.

Moment matrix. Given a sequence $\mathbf{y}=(y_{\alpha})_{\alpha\in\mathbb{N}^{n}}$, let $L_{\mathbf{y}}:\mathbb{R}[\mathbf{x}]\to\mathbb{R}$ be the linear (Riesz) functional

[TABLE]

Given $\mathbf{y}$ and $d\in\mathbb{N}$, the moment matrix associated with $\mathbf{y}$, is the real symmetric $s(d)\times s(d)$ matrix $\mathbf{M}_{d}(\mathbf{y})$ with rows and columns indexed in $\mathbb{N}^{n}_{d}$ and with entries

[TABLE]

Localizing matrix. Given a sequence $\mathbf{y}=(y_{\alpha})_{\alpha\in\mathbb{N}^{n}}$, and a polynomial $g\in\mathbb{R}[\mathbf{x}]$, the localizing moment matrix associated with $\mathbf{y}$ and $g$, is the real symmetric $s(d)\times s(d)$ matrix $\mathbf{M}_{d}(g\,\mathbf{y})$ with rows and columns indexed in $\mathbb{N}^{n}_{d}$ and with entries

[TABLE]

II-B The volume of a compact semi-algebraic set

In this section we recall how to approximate as closely as desired the Lebesgue volume of a compact semi-algebraic set $\mathbf{K}\subset\mathbb{R}^{n}$. It will be the building block of the methodology to approximate the set $\mathbf{X}^{*}_{\epsilon}$ in (I).

Let $\mathbf{X}\subset\mathbb{R}^{n}$ be a box and let $\lambda\in M(\mathbf{X})$ be the Lebesgue measure on $\mathbf{X}$. Let $\mathbf{K}:=\{\mathbf{x}:g_{j}(\mathbf{x})\geq 0,\>j=1,\ldots,m\}$, assumed to be compact. For convenience and with no loss of generality we may and will assume that $g_{1}(\mathbf{x})=M-\|\mathbf{x}\|^{2}$ for some $M>0$.

Theorem II.1** ([5])**

Let $\mathbf{K}\subset\mathbf{X}$ and with nonempty interior. Then

[TABLE]

and $d\phi^{*}=1_{\mathbf{K}}(\mathbf{x})d\lambda$ is the unique optimal solution.

Problem (7) is an infinite-dimensional LP with dual

[TABLE]

where $\mathcal{C}(\mathbf{X})$ is the space of continuous functions on $\mathbf{X}$. That (8) and (9) have same optimal value follows from compactness of $\mathbf{X}$ and Stone-Weierstrass Theorem. Next, as shown in [5], there is no duality gap, i.e., $\rho=\lambda(\mathbf{K})$.

II-C Semidefinite relaxations

Let $d_{j}=\lceil{\rm deg}(g_{j})/2\rceil$, $j=1,\ldots,m$. To approximate $\lambda(\mathbf{K})$ one solves the hierarchy of semidefinite programs, indexed by $d\in\mathbb{N}$:

[TABLE]

Interpretation.

Ideally, the variables $\mathbf{y}=(y_{\alpha})$ (resp. $\mathbf{z}=(z_{\alpha})$) of (10) should be viewed as “moments” of the measure $\phi$ (resp. the measure $\psi:=\lambda-\phi$) in (7) (and so $L_{\mathbf{y}}(1)=\phi(\mathbf{K})$); the constraints $\mathbf{M}_{d}(\mathbf{y}),\mathbf{M}_{d-d_{j}}(g_{j}\,\mathbf{y})\succeq 0$ (resp. $\mathbf{M}_{d}(\mathbf{z})\succeq 0$) are precisely necessary conditions for the above statement to be true111As $\mathbf{X}=[-1,1]^{n}=\{\mathbf{x}:1-x_{j}^{2}\geq 0,\>j=1,\ldots,n\}$, in principle one should also impose $\mathbf{M}_{d-1}((1-x_{j}^{2})\,\mathbf{z})\succeq 0$, $j=1,\ldots,n$, for all $d$, to ensure that $\psi$ is supported on $\mathbf{X}$. However as $\mathbf{M}_{d}(\mathbf{z})\preceq\mathbf{M}_{d}(\lambda)$ for all $d$, in the limit as $d\to\infty$, one has $\psi\leq\mu$ and so ${\rm support}(\psi)\subset\mathbf{X}$. (and which become sufficient as $d\to\infty$).

The sequence $(\rho_{d})_{d\in\mathbb{N}}$ is monotone non increasing and $\rho_{d}\to\lambda(\mathbf{K})$ as $d\to\infty$. However the convergence is rather slow and in [5] the authors have proposed to replace the criterion $L_{\mathbf{y}}(1)$ by $L_{\mathbf{y}}(h)$ where $h\in\mathbb{R}[\mathbf{x}]$ is a polynomial that is nonnegative on $\mathbf{K}$ and vanishes on the boundary of $\mathbf{K}$. If one denotes by $\mathbf{y}^{d}$ an optimal solution of (10) then $\rho_{d}\to\int_{\mathbf{K}}hd\lambda$ and $y^{d}_{0}\to\lambda(\mathbf{K})$ as $d\to\infty$. The convergence $y^{d}_{0}\to\lambda(\mathbf{K})$ is much faster but is not monotone anymore, which can be annoying because we do not obtain a decreasing sequence of upper bounds on $\lambda(\mathbf{K})$ as was the case with (10). For more details the interested reader is referred to [5].

II-D Stokes can help

This is why we propose another technique to accelerate the convergence $\rho_{d}\to\lambda(\mathbf{K})$ in (10) while maintaining its monotonicity. So let $h\in\mathbb{R}[\mathbf{x}]$ be such that $h(\mathbf{x})=0$ for all $\mathbf{x}\in\partial\mathbf{K}$ (but $h$ is not required to be nonnegative on $\mathbf{K}$). Then by Stokes’ theorem (with vector field $X=e_{i}\in\mathbb{R}^{n}$, $e_{ij}=\delta_{i=j}$, $i,j=1,\ldots,n$), for each $\alpha\in\mathbb{N}^{n}$:

[TABLE]

and so the optimal solution $\phi^{*}$ of Theorem II.1 must satisfy

[TABLE]

Therefore in (10) we may impose the additional moment constraints

[TABLE]

To appreciate the impact of such additional constraints on the convergence $\rho_{d}\to\lambda(\mathbf{K})$, consider the simple example with $n=2$ and $\mathbf{X}=[-a,a]^{2}$, let $\mathbf{K}:=\{\mathbf{x}:\|\mathbf{x}\|^{2}\leq 1\}$ so that $\lambda(\mathbf{K})=\pi$. For different values of $a$ and $d=3,4$, results are displayed in Table I.

Remark II.2

Theorem II.1 is valid for any measure $\mu\in M(\mathbf{X})$ and not only the Lebesgue measure $\lambda$. On the other hand, additional Stokes constraints similar to (11) (but now with vector field $X=\mathbf{x}$ and polynomial $\theta_{\alpha}$ below) are valid provided that $d\mu=f\,d\lambda$ or $d\mu=\exp(f)\,d\lambda$ for some homogeneous polynomial $f\in\mathbb{R}[\mathbf{x}]$. Then with $d_{f}={\rm deg}(f)$,

[TABLE]

III Main Result

After the preliminary results of §II, we are now in position to state our main result. Let $\mu$ be the distribution of the noise parameter $\mathbf{\omega}\in\mathbf{\Omega}$, and let $\lambda$ be the Lebesgue measure on $\mathbf{X}$. The notation $\lambda\otimes\mu$ denotes the product measure on $\mathbf{X}\times\mathbf{\Omega}$, that is,

[TABLE]

With $\mathbf{K}\subset\mathbf{X}\times\mathbf{\Omega}$ as in (1), and for every $\mathbf{x}\in\mathbf{X}$, let $\mathbf{K}_{\mathbf{x}}$ be as in (2) (possibly empty). Consider the infinite dimensional LP:

[TABLE]

Theorem III.1

The unique optimal solution of (12) is

[TABLE]

and the optimal value $\rho$ of (12) satisfies

[TABLE]

Proof:

That $\rho=\displaystyle\int_{\mathbf{X}\times\mathbf{\Omega}}1_{\mathbf{K}}((\mathbf{x},\mathbf{\omega}))\,\lambda\otimes\mu(d(\mathbf{x},\mathbf{\omega}))$ follows from Theorem II.1 (with $\lambda\otimes\mu$ instead of $\lambda$ in Theorem II.1). By Fubini-Tonelli’s Theorem (see Ash [1][Theorem 2.6.6, p. 103])

[TABLE]

[TABLE]

∎

Semidefinite relaxations

Let $d_{j}=\lceil{\rm deg}(g_{j})/2\rceil$ for all $j$. As we did for (10) in §II, let $\mathbf{y}=(y_{\alpha,\beta})$ and $\mathbf{z}=(y_{\alpha,\beta})$, $(\alpha,\beta)\in\mathbb{N}^{n+p}$, and relax (12) to the following hierarchy of semidefinite programs, indexed by $d\in\mathbb{N}$:

[TABLE]

and of course $\rho_{d}\geq\rho_{d+1}\geq\rho$ for all $d$. The dual of (14) is the semidefinite program:

[TABLE]

Again as $\mathbf{K}$ is compact, for technical reasons (but with no loss of generality) we may and will assume that in the definition (1) of $\mathbf{K}$, $g_{1}(\mathbf{x})=M-\|\mathbf{x}\|^{2}$ for some $M>0$.

Theorem III.2

Let $\mathbf{K}$ and $(\mathbf{X}\times\mathbf{\Omega})\setminus\mathbf{K}$ be with nonempty interior. There is no duality gap between (14) and its dual (15), i.e., $\rho_{d}=\rho^{*}_{d}$ for all $d$. In addition (15) has an optimal solution $p^{*}_{d}\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]_{2d}$ such that

[TABLE]

Define $h^{*}_{d}\in\mathbb{R}[\mathbf{x}]_{2d}$ to be:

[TABLE]

Then $h^{*}_{d}(\mathbf{x})\geq\mu(\mathbf{K}_{\mathbf{x}})$ for all $\mathbf{x}\in\mathbf{X}$ and

[TABLE]

as $d\to\infty$.

Proof:

That $\rho_{d}=\rho^{*}_{d}$ is because Slater’s condition holds for (14). Indeed let $\mathbf{y}^{*}$ be the moments of $\phi^{*}$ in Theorem III.1 and $\mathbf{z}^{*}$ be the moments of $\lambda\otimes\mu-\phi^{*}$ (on $(\mathbf{X}\times\mathbf{\Omega})\setminus\mathbf{K}$). Then as $\mathbf{K}$ has nonempty interior, $\mathbf{M}_{d}(\mathbf{y}^{*})\succ 0$ and $\mathbf{M}_{d}(g_{j}\,\mathbf{y}^{*})\succ 0$ for all $d$. Similarly as $(\mathbf{X}\times\mathbf{\Omega})\setminus\mathbf{K}$ has nonempty interior, $\mathbf{M}_{d}(\mathbf{z}^{*}))\succ 0$. Moreover since the optimal value $\rho_{d}$ is finite for all $d$ this implies that (15) has an optimal solution $p^{*}_{d}\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]_{2d}$. Therefore:

[TABLE]

where $h^{*}_{d}(\mathbf{x})\geq\mu(\mathbf{K}_{\mathbf{x}})$ follows from $p^{*}_{d}\geq 1$ on $\mathbf{K}$. Finally the convergence $\lim_{d\to\infty}\rho_{d}=\rho$ follows from Theorem II.1. ∎

Then as $h^{*}(\mathbf{x})\geq\mu(\mathbf{K}_{\mathbf{x}})$ on $\mathbf{X}$, the sets $\mathbf{X}^{d}_{\epsilon}=\{\mathbf{x}\in\mathbf{X}:h^{*}_{d}(\mathbf{x})>1-\epsilon\}$, $d\in\mathbb{N}$, form a sequence of outer approximations of the set $\mathbf{X}^{*}_{\epsilon}$. In fact more can be said.

Corollary III.3

Let $h^{*}_{d}\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]_{2d}$ be as in Theorem III.2. Then the function $\mathbf{x}\mapsto\psi^{*}_{d}(\mathbf{x}):=h^{*}_{d}(\mathbf{x})-\mu(\mathbf{K}_{\mathbf{x}})$ is nonnegative on $\mathbf{X}$ and converges to [math] in $L_{1}(\mathbf{X},\lambda)$. In particular $\psi^{*}_{d}\to 0$ in $\lambda$-measure222A sequence of functions $h,(h_{n})_{n\in\mathbb{N}}$ on a measure space $(\mathbf{X},\mathcal{B},\lambda)$ converges in measure to $h$ if for all $\epsilon>0$, $\lambda(\{\mathbf{x}:|h_{n}(\mathbf{x})-h(\mathbf{x})|\geq\epsilon\})\to 0$ as $n\to\infty$. The sequence converges almost-uniformly to $h$ if to every $\epsilon>0$ there exists a set $B_{\epsilon}\subset\mathbf{X}$ such that $\lambda(B_{\epsilon})<\epsilon$ and $h_{n}\to h$ uniformly in $\mathbf{X}\setminus B_{\epsilon}$., and $\lambda$-almost uniformly for some subsequence $(\psi^{*}_{d_{k}})_{k\in\mathbb{N}}$.

Proof:

As $\rho_{d}\to\rho$ as $d\to\infty$,

[TABLE]

whence the convergence to [math] in $L_{1}(\mathbf{X},\lambda)$. Then convergence $\psi^{*}_{d}\to 0$ in $\lambda$-measure, and $\lambda$-almost sure convergence for a subsequence follow from standard results from Real Analysis. See e.g. Ash [1, Theorem 2.5.1]. ∎

As we next see, the convergence $h^{*}_{d}(\mathbf{x})\to\mu(\mathbf{K}_{\mathbf{x}})$ in $\lambda$-measure established in Corollary III.3 will be useful to obtain strong asymptotic guarantees.

III-A Strong asymptotic guarantees

We here investigate asymptotic properties of the sequence of sets $(\mathbf{X}^{d}_{\epsilon})_{d\in\mathbb{N}}$, as $d\to\infty$.

Corollary III.4

With $\mathbf{X}^{*}_{\epsilon}$ as in (I), let $\mathbf{X}^{d}_{\epsilon}:=\{\mathbf{x}\in\mathbf{X}:h^{*}_{d}(\mathbf{x})\geq 1-\epsilon\}$ where $h^{*}_{d}$ is as in Theorem III.2, $d\in\mathbb{N}$. Then:

[TABLE]

Proof:

Observe that

[TABLE]

and therefore

[TABLE]

Next, for each $\ell=1,\ldots$, write

[TABLE]

By the convergence $h^{*}_{d}\to\mu(\mathbf{K}_{\mathbf{x}})$ in $\lambda$-measure as $d\to\infty$, $\lim_{d\to\infty}\lambda(R_{\ell}\cap\mathbf{X}^{d}_{\epsilon})=0$ and so

[TABLE]

This implies

[TABLE]

which in turn yields the desired result (16). ∎

III-B Inner Approximations

In the previous section we have provided a converging sequence $(\mathbf{X}^{d}_{\epsilon})_{d\in\mathbb{N}}$ of outer approximations of $\mathbf{X}^{*}_{\epsilon}$. Clearly, letting $\mathbf{K}^{c}:=(\mathbf{X}\times\mathbf{\Omega})\setminus\mathbf{K}$, the same methodology now applied to a chance constraint of the form

[TABLE]

would provide a converging sequence of inner approximations of the set $\mathbf{X}^{*}_{\epsilon}:=\{\mathbf{x}\in\mathbf{X}:{\rm Prob}((\mathbf{x},\mathbf{\omega})\,\in\mathbf{K})\,>\,1-\epsilon\}$. To do so, (i) write $\mathbf{K}^{c}$ as a finite union $\bigcup_{i}\mathbf{K}^{c}_{i}$ of basic semi-algebraic sets $\mathbf{K}^{c}_{i}$ (whose $\mu\otimes\lambda$ measure of their overlaps is zero), (iii) apply the above methodology to each $\mathbf{K}^{c}_{i}$, and then (ii) sum-up the results.

III-C Accelerating convergence

As we already have seen in Section II-D for the semidefinite program (10), as $d\to\infty$ the convergence $\rho_{d}\to\rho$ of the optimal value of (14) can also be slow due to the Gibb’s phenomenon333The Gibbs’ phenomenon appears at a jump discontinuity when one approximates a piecewise $C^{1}$ function with a continuous function, e.g. by its Fourier series. that appears in the dual (15) when approximating the indicator function $\mathbf{x}\mapsto 1_{\mathbf{K}}(\mathbf{x})$ by a polynomial.

So assume that $\mu$ is the Lebesgue measure on $\mathbf{\Omega}$ where for instance $\mathbf{\Omega}=[-1,1]^{p}$, scaled to be a probability measure (but the same idea works if $d\mu=h(\mathbf{\omega})d\mathbf{\omega}$, or if $d\mu=\exp(h(\mathbf{\omega}))d\mathbf{\omega}$ for s ome homogeneous polynomial $h$). Then again we propose to include additional constraints on the moments $\mathbf{y}$ and $\mathbf{z}$ in (14) coming from additional properties of the optimal solution $\phi^{*}$ and $\psi^{*}=\lambda\otimes\mu-\phi^{*}$ of (12). Again these additional properties are coming from Stokes’ formula but now for integrals on $\mathbf{K}_{\mathbf{x}}$ (resp. $\mathbf{\Omega}\setminus\mathbf{K}_{\mathbf{x}}$), then integrated over $\mathbf{X}$.

Let $f_{1},f_{2}\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]$ be the polynomials $(\mathbf{x},\mathbf{\omega})\mapsto f_{1}(\mathbf{x},\mathbf{\omega}):=\prod_{j}g_{j}(\mathbf{x},\mathbf{\omega})$ and $(\mathbf{x},\mathbf{\omega})\mapsto f_{2}(\mathbf{x},\mathbf{\omega}):=f_{1}(\mathbf{x},\mathbf{\omega})\prod_{j}(1-\mathbf{\omega}_{j}^{2})$ of respective degree $d_{1}$, $d_{2}$. For each fixed $\mathbf{x}\in\mathbf{X}$, the polynomial $\mathbf{\omega}\mapsto f_{1}(\mathbf{x},\mathbf{\omega})$ (resp. $\mathbf{\omega}\mapsto f_{2}(\mathbf{x},\mathbf{\omega})$) vanishes on the boundary $\partial\mathbf{K}_{\mathbf{x}}$ of $\mathbf{K}_{\mathbf{x}}$ (resp. $\partial\mathbf{K}_{\mathbf{x}}^{c}$ of $\mathbf{\Omega}\setminus\partial\mathbf{K}_{\mathbf{x}}$). Therefore for each $\beta\in\mathbb{N}^{p}$, Stokes’ Theorem (applied with vector fields $e_{j}\in\mathbb{R}^{p}$ (where $e_{jk}=(\delta_{j=k})$), $k,j=1,\ldots,p$), states:

[TABLE]

[TABLE]

So let $\theta^{k}_{j\beta}\in\mathbb{R}[\mathbf{x},\mathbf{\omega}]$ of degree $d_{k}+|\beta|-1$, $k=1,2$, be:

[TABLE]

for all $\beta\in\mathbb{N}^{p}$, $j=1,\ldots,p$. Then for each $(\alpha,\beta)\in\mathbb{N}^{n+p}$:

[TABLE]

[TABLE]

Equivalently, in view of what are $\phi^{*},\psi^{*}$ in Theorem III.1,

[TABLE]

Therefore in (14) we may include the additional moments constraints $L_{\mathbf{y}}(\mathbf{x}^{\alpha}\,\theta^{1}_{j\beta}(\mathbf{x},\mathbf{\omega}))=0$, and $L_{\mathbf{z}}(\mathbf{x}^{\alpha}\,\theta^{2}_{j\beta}(\mathbf{x},\mathbf{\omega}))=0$, for all $(\alpha,\beta)\in\mathbb{N}^{n+p}$ such that $|\alpha+\beta|\leq 2d+1-d_{1}$ and $|\alpha+\beta|\leq 2d+1-d_{2}$ respectively.

III-D The non-compact case

In some applications the noise $\mathbf{\omega}$ is assumed to follow a normal distribution $\mu$ on $\mathbf{\Omega}=\mathbb{R}^{p}$. Therefore $\mathbf{\Omega}$ is not compact anymore and the machinery used in [5] cannot be applied directly. However the normal distribution satisfies the important Carleman’s property. That is, let $L_{\mathbf{y}}$ be the Riesz functional associated with $\mu$, i.e., $L_{\mathbf{y}}(f)=\int fd\mu$ for all $f\in\mathbb{R}[\mathbf{\omega}]$. Then

[TABLE]

In particular $\mu$ is moment determinate, that is, $\mu$ is completely determined by its moments. These two properties have been used extensively in e.g. Lasserre [11] and also in [12], precisely to show that with $\mathbf{K}\subset\mathbf{\Omega}$ not necessarily compact, one may still approximate its Gaussian measure $\mu(\mathbf{K})$ as closely as desired. Again one solves (7) via the same hierarchy of semidefinite relaxations (10) (but now with $\mu$ instead of $\lambda$). For more details the interested reader is referred to Lasserre [11, 12].

In view of the above (technical) remarks, one may then extend the machinery described in §III to the case where $\mathbf{\Omega}=\mathbb{R}^{p}$, $\mu$ is the Gaussian measure, and $\mathbf{K}_{\mathbf{x}}$ ($\mathbf{x}\in\mathbf{X}$) is not necessarily compact. A version of Stokes’ Theorem for non compact sets is even described in [12] to accelerate the convergence of the semidefinite relaxations (10) (with $\mu$ instead of $\lambda$). It can be used to accelerate the convergence of the semidefinite relaxations (14), exactly as we do in §III-C for the compact case.

III-E Numerical examples

For illustration purposes we have considered simple small dimensional examples for which the function $\mathbf{x}\mapsto\mu(\mathbf{K}_{\mathbf{x}})$ has a closed form expression, so that we can compare the set $\mathbf{X}^{*}_{\epsilon}$ with its approximations $\mathbf{X}^{d}_{\epsilon}$, $d\in\mathbb{N}$, obtained in Corollary III.4 (with and without using Stokes constraints).

Example III.5

$\mathbf{X}=[-1,1]$, $\mathbf{\Omega}=[0,1]$ and $\mathbf{K}=\{(x,\omega):1-x^{2}/0.81-\omega^{2}/1.44\geq 0\}$. $\lambda$ and $\mu$ are the Lebesgue measure. In this case $\mathbf{X}^{*}_{\epsilon}=[a_{\epsilon},b_{\epsilon}],\mathbf{X}^{d}_{\epsilon}=[a^{d}_{\epsilon},b^{d}_{\epsilon}]\subset[-1,1]$. In Table II we display the relative error $(b^{d}_{\epsilon}-b_{\epsilon}+a_{\epsilon}-a^{d}_{\epsilon})/(b^{d}_{\epsilon}-a^{d}_{\epsilon})$ for different values of $\epsilon$ and $d$, with and without Stokes constraints. The results indicate that adding Stokes constraints help a lot. With relatively few moments $2d\leq 16$ one obtains good approximations.

Example III.6

$\mathbf{X}=[-1,1]$, $\mathbf{\Omega}=[0,1]$ and $\mathbf{K}=\{(x,\omega):1-x^{2}-\frac{\omega^{2}}{2}\geq 0;\>\frac{x^{2}}{2}+\omega^{2}-\frac{1}{4}\geq 0\}$. $\lambda$ and $\mu$ are the Lebesgue measure on $\mathbf{X}$ and $\mathbf{\Omega}$ respectively. In this case, when $\epsilon<0.4$, the set $\mathbf{X}^{*}_{\epsilon}$ is the union of two disconnected intervals, hence more difficult to approximate. As in Example III.5, in Table III we display the relative error for different values of $\epsilon$ and $d$, with and without Stokes constraints. Again, the results indicate that adding the Stokes constraints help a lot. With relatively few moments $2d\leq 20$ one obtains good approximations. For instance, $\mathbf{X}^{*}_{0.1}=[-0.7714,-0.3082]\cup[0.3082,0.7714]$ while one obtains $\mathbf{X}^{*}_{0.1}\subset\mathbf{X}^{10}_{0.1}=[-0.7985,-0.26]\cup[0.26,0.7985]$ with Stokes and the larger set $\mathbf{X}^{10}_{0.1}=[-0.8673,-0.1881]\cup[0.1881,0.8673]$ without Stokes constraints.

Example III.7

$\mathbf{X}=[-1,1]^{2}$, $\mathbf{\Omega}=[0,1]$ and $\mathbf{K}=\{(x,\omega):1-x_{1}^{2}-x_{2}^{2}-\omega^{2}\geq 0\}$. $\lambda$ and $\mu$ are the Lebesgue measure on $\mathbf{X}$ and $\mathbf{\Omega}$ respectively. For this two-dimensional example (in $\mathbf{x}$) we have plotted the boundary of the set $\mathbf{X}^{*}_{\epsilon}$ (inner approximate circle, solid line in black). The curve in the middle (red dashed line) (resp. outer, blued dashed line) is the boundary of the approximation $\mathbf{X}^{d}_{\epsilon}$ computed with Stokes constraints (resp. without Stokes constraints). For $\epsilon=0.01$ and $d=10$ the results are plotted in Fig. 1 and in Fig. 2 for $\epsilon=0.05$ and $d=10$.

IV CONCLUSION

We have presented a systematic numerical scheme to provide an sequence of outer and inner approximations $(\mathbf{X}^{d}_{\epsilon})$ of the feasible set $\mathbf{X}^{*}_{\epsilon}$ associated with chance constraints of the form ${\rm Prob}((\mathbf{x},\mathbf{\omega})\in\mathbf{K})>1-\epsilon$. Each outer and inner approximation $\mathbf{X}^{d}_{\epsilon}$ is the [math]-super level set of some polynomial whose coefficients are computed via solving a certain semidefinite program. As $d$ increases $\lambda(\mathbf{X}^{d}_{\epsilon}\setminus\mathbf{X}^{*}_{\epsilon})\to 0$, a nice and highly desirable asymptotic property. Of course this methodology is computationally expensive and in its present form limited to problems of small size. But we hope it can pave the way to define efficient heuristics. Also checking whether this methodology can accommodate potential sparsity patterns present in larger size problems, is a topic of further investigation.