# Soft-constrained interval predictor models and epistemic reliability intervals: a new tool for uncertainty quantification with limited experimental data

Roberto Rocchetta, Qi Gao, Milan Petkovic

► To cite this version:

Roberto Rocchetta, Qi Gao, Milan Petkovic. Soft-constrained interval predictor models and epistemic reliability intervals: a new tool for uncertainty quantification with limited experimental data. Mechanical Systems and Signal Processing, 2021, 161, pp.107973. 10.1016/j.ymssp.2021.107973. hal-04457309

HAL Id: hal-04457309

https://minesparis-psl.hal.science/hal-04457309v1

Submitted on 14 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Soft-constrained interval predictor models and epistemic reliability intervals: a
new tool for uncertainty quantification with limited experimental data

Roberto Rocchettaª, Qi Gaoᵇ, Milan Petkovicª,ᵇ

ªSecurity group, Mathematics and Computer Science,
Technical University of Eindhoven, The Netherlands

<sup>b</sup>Department of Data Science, Philips Research,
Eindhoven, The Netherlands

**Abstract**

Interval Predictor Models (IPMs) offer a non-probabilistic, interval-valued, characterization of the uncertainty affecting random data generating processes. IPMs are constructed directly from data, with no assumptions on the distributions of the uncertain factors driving the process, and are therefore exempt from the subjectivity induced by such a practice. The reliability of an IPM defines the probability of correct predictions for future samples and, in practice, its true value is always unknown due to finite samples sizes and limited understanding of the process. This paper proposes an overview of scenario optimization programs for the identification of IPMs. Traditional IPM identification methods are compared with a new scheme which softens the scenario constraints and exploits a trade-off between reliability and accuracy. The new methods allows prescribing predictors that achieve higher accuracy for a quantifiable reduction in the reliability. Scenario optimization theory is the mathematical tool used to prescribe formal epistemic bounds on the predictors reliability. A review of relevant theorems and bounds is proposed in this work. Scenario-based reliability bounds hold distribution-free, non asymptotically, and quantify the uncertainty affecting the model’s ability to correctly predict future data. The applicability of the new approach is tested on three examples: i) on the modelling of a trigonometric function affected by a noise term, ii) on the identification of a black-box system-controller dynamic response model and, iii) on the modelling of the vibration response of a car suspension arm crossed by a crack of unknown length. The points of strength and limitations of the new IPM are discussed based on the accuracy, computational cost, and width of the resulting epistemic bounds.

***Keywords:*** Interval predictor models, Scenario theory, Epistemic uncertainty, Reliability bounds, Convex Optimization, Frequency response function

# 1. Introduction

An Interval Predictor Model is a rule that assigns to explanatory inputs ranges of plausible outcomes of a dependent variable [1]. The uncertainty in the process is quantified via non-probabilistic, interval-valued, predictions [2]. Uncertainties can be modelled probabilistically, non-probabilistically or by a mixture of these two approaches [3, 4, 5]. Probabilistic approaches often require unwarranted assumptions to model epistemic uncertainty and this can lead to underestimation of the true level of uncertainty if a lack of data affects the study [6]. Most importantly, a decision taken based on a wrong characterization may be biased, misleading, and potentially harmful. If a large epistemic uncertainty affects our study, the use of a mixed approach or non-probabilistic approach is advisable [7]. Evidence theory [8, 9], Possibility theory [10], Credal sets, Fuzzy sets and Ambiguity sets theory [11, 12, 13, 14], are some of the most commonly used paradigms for this [15].

In contrast to traditional methods, Interval Predictor Models do not require distributional assumptions to describe the data and are exempt from the subjectivity induced by such a practice. IPMs are constructed directly from data by fitting two functions to envelop a realization of a random Data Generating Mechanism (DGM). The parameters of these functions are optimized by imposing sample-constraints and enforcing a tightness condition on the width. This makes an IPM a data-driven and distribution-free model [16, 17]. A robust IPM must give both accurate and reliable predictions. The accuracy of an IPM is the average width of the bounds whilst the IPM reliability is the probability of a new sample correctly falling within an interval prediction. Accuracy and reliability are competitive

Email addresses: r.rocchetta@tue.nl (Roberto Rocchetta), q.gao@philips.com (Qi Gao), milan.petkovic@philips.com (Milan Petkovic)
objectives since wider intervals are more likely to contain new data but lead to less accurate prediction. Unfortunately, reliability is only known approximatively via samples-based estimations and, in practice, its exact value is always unavailable due to finite sample sizes. Scenario optimization theory offers a mathematical framework to quantify the uncertainty in the reliability by prescribing formal epistemic bounds on its real (unknown) value. The width of the bounds quantifies the effect of epistemic uncertainty.

Similarly to statistical learning theory [18, 19] and compression learning [20, 21], Scenario theory equips the solution of learning problems with a probability of error certificate. Within scenario theory, a learning problem is a *scenario optimization program*, i.e., a program where a finite number of samples (*scenarios*) enforce deterministic constraints in the optimization. Scenario theory has been extensively studied for hard-constrained convex programs, [22, 23, 24, 25, 26], and was recently extended to non-convex and soft-constrained cases, [27, 28, 29, 30]. For the solution of a convex problem, Campi and Garatti [24] showed that the error probability distribution is bounded by a Beta distribution whose parameters only depend on the number of available samples and the dimension of the optimization problem. This bound results tight (exact) for the prototype class of fully-supported problems [24], i.e., problems for which the decision complexity is equal to the dimension of the design space with probability one. More in general, Carè et al. [31] proven that the distribution of empirical reliability levels follows an ordered Dirichlet distribution whose marginals are Beta distributions.

Scenario theory has been used to solve several problems in different disciplines. As example, a scenario approach was used to tackle uncertainty characterization problems [32], a water flooding optimization [33], an antenna arrays design [34], optimal control and reliability-based design-optimization problems [35, 36, 37]. More recently, the generalization of Support Vector Machine models [38], Game-theoretic policies [39] and other machine learning tools [40, 41, 42, 43, 44] were also assessed via Scenario theory. Scenario optimization theory is often used to quantify the uncertainty affecting the reliability of IPMs. As example, IPMs with guaranteed reliability were used for history matching problems [45], space radiation shielding analysis [46], damaged beams dynamic response modeling [47], design of adaptive controller for Aerial Vehicle's [48], robust predictive control [49], robust uncertainty propagation [50], and for the detection and removal outliers [2]. In [51] a random predictor model was introduced and compared with IPMs. The interested reader is reminded to, e.g., [36, 52, 53] for further reading on efficient algorithmic implementations of Scenario theory and dedicated software. For a detailed discussion on the connection with statistical learning theory and compression learning theory see, for instance, [18, 20, 21, 29].

This paper presents an overview of scenario optimization programs and test their use for prescribing interval predictors. A survey of scenario theoretic epistemic bounds is also proposed, relevant theorems are presented and discussed based on the numerical results. A new IPM optimization method, named soft-constrained IPM program, is originally introduced in this work. The method based on the constraints relaxation scheme studied by [54] allows trading model accuracy for reliability. The accuracy and reliability of the new method are tested on three realistic problems: i) on the modelling of an unknown trigonometric function affected by noise [55]; ii) on the identification of a model for the dynamic response of a black-box system-controller affected by uncertainty [56]; iii) on the modelling of the frequency response function of a car suspension arm crossed by a crack of unknown length. Case study iii) is new and is introduced in this work based on the simulated data presented in [57].

The main contributions of this work can be summarized as follows:

1. We present an overview of scenario programs for the identification of Interval Predictor Models.

2. We present scenario-based epistemic bounds and compare them numerically on easily reproducible examples.

3. We propose a new soft-constrained IPM identification methods. The method can be used to trade reliability for accuracy in the predictions.

4. We promote the applicability of IPMs to model the response of two realistic dynamical systems affected by uncertainty. In particular, the modelling of the frequency response of a car suspension arm crossed by a crack of unknown length is originally introduced in this work.

The rest of this paper is organized as follows: Section 2 introduces optimization programs for interval predictor models and discusses model accuracy and reliability. Scenario optimization programs and reliability bounds are reviewed in Section 3. Section 4 introduces a soft-constrained optimization program for IPMs. Section 5 presents the case studies, the numerical results, and discussions. Section 6 closes the paper.

## 2. Interval Predictor Model

A stochastic data generating mechanism is defined by a probability space $(\Delta, \mathcal{F}, \mathbb{P})$, where $\Delta$ is the event space equipped with a $\sigma$-algebra $\mathcal{F}$ and a probability measure $\mathbb{P}$. The DGM generates samples $(x(t), y(t))$ from an unknown random process $y(t) = f(x(t)) + \mathcal{E}$, where $\mathcal{E}$ is a random noise term and a pair $\delta(t) = (x(t), y(t)) \in \Delta$ is composed
by a vector of explanatory variables $x(t) \in \mathcal{X} \subseteq \mathbb{R}^{n_x}$ and a vector of dependent variables $y(t) \in \mathcal{Y} \subseteq \mathbb{R}^{n_y}$ such that $\Delta = \mathcal{X} \times \mathcal{Y}$. A DGM represents, for instance, the processes of collecting measurements from a dynamical system affected by uncertainty or gathering numerical experiments from a stochastic simulator (a physics-based digital twin of the real system). The function $f(x(t))$ and noise distribution are unknown in practice, and only a set of $N$ realizations $\mathcal{D}_N = \{x_i, y_i\}_{i=1}^N$ available for the characterization of the process. The input-output samples $\delta_i = (x_i, y_i) \in \mathcal{D}_N$ will be referred to as scenarios.

Traditional regression analysis can be used to identify a deterministic model which is optimal in the least square error sense, that is, a function $\hat{f}(x) : \mathbb{R}^{n_x} \to \mathbb{R}^{n_y}$, approximating the true $f(x)$. Although $\hat{f}(x)$ can capture well the overall behavior of the unknown $f(x)$, the uncertainty affecting the process and the data limits the predictive ability $\hat{f}(x)$. To overcome this limitation, a non-deterministic regression model can be considered which accounts for the relevant uncertainties by replacing the one-input to a one-output function $\hat{f}(x)$, with a map from one-input to many-outputs. Non-deterministic regression models can be probabilistic or non-probabilistic depending on how the uncertainty is modelled, traditional examples are Gaussian processes [58, 59, 60, 61, 62], Wiener and Poisson processes [63], Interval-valued B-splines [64], support vector regression models [65, 66], neural networks [67, 68], random fields [69], Fuzzy regression models [70], interval response surfaces [71], or Random predictive models [51, 72].

Interval Predictive Model offers an interval-valued characterization of the uncertainty affecting a stochastic data generating process [46] by defining a rule $\mathcal{I}_N(x; \theta)$, which assigns to a vector explanatory variable $x$ a lower and upper bound $\mathcal{I}_N(x; \theta) = [f_l(x; \theta), f_u(x; \theta)]$. A vector of fitting coefficients $\theta$ defines two functions $f_l(x; \theta)$ and $f_u(x; \theta)$ bounding the dependent quantity of interest $y$. The resulting interval $\mathcal{I}_N(x; \theta)$, gives a non-probabilistic characterization of the uncertainty because a random realization $(x, y)$ is considered to be unknown except that it belongs to a given set $y \in \mathcal{I}_N(x; \theta)$. The reminder of this work will focus on scalar input-output models to ease the presentation of the approach and visualization of the results, i.e., $n_x = n_y = 1$ and $n_\delta = 2$. This has not to be considered as a limiting assumption. In fact, the dimension $n_\delta$ has no impact on the applicability of the approach to higher dimensions nor (in principle) on the computational cost of IPM optimization methods which will be presented next.

## 2.1. Stationary and independent process

For simplicity sake, we assume the process $\delta_i = (x_i, y_i)$, $i = 1, 2, ...$ to be a stationary and independent sequence. If the subscripts $i$ defines a time index $t$, the marginal distribution function of a stationary process at time $t$ is the same at any $t$:

$$F_{\delta}(\delta_{1+\tau}, \dots, \delta_{l+\tau}, \dots) = F_{\delta}(\delta_1, \dots, \delta_l, \dots) \quad \forall \tau, t \in \mathbb{R}.$$

In other words, a stochastic process is stationary if the distribution of $\delta(t)$ is invariant in time, that is, the distribution of the random process is the same for any $t$. Differently, the independence assumption on the process $\delta_i$, $t = 1, 2, ...$ rules out the presence of inter-time correlations:

$$F_{\delta_1, \dots, \delta_i}(\delta_1, \dots, \delta_i) = F_{\delta_1}(\delta_1) \cdots F_{\delta_i}(\delta_i) \quad \forall \delta_i.$$

Independence in $x(t)$ is an assumption which is generally met in many identification problems such as in economic, pattern recognition, learning theory where no dynamics is present [16]. On the other hand, it may be of interest to generalize the IPM framework to a dependent context. This is still unavailable in a general form and, thus, the application of IPM to data dependent sequences in full generality require further research. Nonetheless, the theory presented in this work is approximately applicable to correlated processes and given that the correlation pattern is estimated and compensated via a deconvolution procedure. For instance, [17] presents a discussion on the applicability of IPM to M-dependent processes. An M-dependent process defines situations where the dependence between points of a stochastic process can be considered negligible after M time steps.

## 2.2. Basis functions

The functions $f_u$ and $f_l$ are represented as a linear combination of fitting coefficients and basis function:

$$f_u(x; \theta) = \sum_{j=0}^{n_u} u_j \psi_j(x), \quad f_l(x; \theta) = \sum_{j=0}^{n_l} l_j \psi_j(x), \qquad (1)$$

where $\theta = [u_0, \dots, u_{n_u}, l_0, \dots, l_{n_l}] \in \mathbb{R}^{n_\theta}$ is the vector of fitting coefficients, $\psi_j(x)$ are basis functions, and $n_u$ and $n_l$ are the number of coefficients for $f_u$ and $f_l$, respectively. The total number of parameters defining an IPM is $n_\theta = n_l + n_u$. Table 1 presents examples of basis function $\psi_j(x)$ which can be used to construct an IPM. Different basis, such as Fuzzy [73] or Wavelets basis [74], can be considered to improve the effectiveness of the model. The problem of selecting a good $\psi_j(x)$ is data and process dependent and in the remainder of this work $\psi_j(x)$ will be selected on an empirical basis.
## 2.3. Accuracy of an IPM

An IPM optimization program seeks an optimal vector $\theta^*$ of fitting coefficients such that the predictor $\mathcal{I}_N(x; \theta^*)$ maximizes the predictive accuracy and reliability of the predictions. The model $\mathcal{I}_N(x; \theta)$ can be optimized to fully enclose the samples in $\mathcal{D}_N$. Clearly, the condition $y_i \in \mathcal{I}_N(x_i; \theta)$ for all $i = 1, \dots, N$ is not a sufficient condition for the predictor $\mathcal{I}_N$ to be accurate. To maximize its accuracy, a tightness condition on $\mathcal{I}_N(x; \theta)$ must be enforced. To this end, the following area metric is minimized:

$$ \mathcal{A}(\theta) = \int_{\mathcal{X}} f_u(x; \theta) - f_l(x; \theta)dx. \qquad (2) $$

The integral in (2) can be approximated via standard numerical integration rules, e.g., midpoint rule or trapezoidal rule. For instance, the area can be obtained by selecting an arbitrary number of points $x_0, x_1, \dots, x_m \in \mathcal{X}$ in descending order and equally spaced by a distance $\Delta_x = x_{j+1} - x_j$. Then, the points are used for the estimation of the integral as follows,

$$ \mathcal{A}(\theta) \approx \Delta_x \sum_{j=1}^{m} \frac{f_u(x_{j-1}; \theta) - f_u(x_j; \theta)}{2} - \Delta_x \sum_{j=1}^{m} \frac{f_l(x_{j-1}; \theta) - f_l(x_j; \theta)}{2}. $$

A large number of integration points $m$ leads to lower approximation errors but also increases the computational cost of the estimation of the integral.

Alternatively, an analytical derivatoin of the integral may be used to speed up the calculations. If the integration domain $\mathcal{X}$ is a one dimensional integral $[\underline{x}, \bar{x}]$ and polynomial basis are selected, the area can be expressed analytically as follows:

$$ \mathcal{A}(\theta) = \sum_{j=0}^{n_u} u_j \left[ \frac{x^{j+1}}{j} \right]_{\underline{x}}^{\bar{x}} - \sum_{j=0}^{n_l} l_j \left[ \frac{x^{j+1}}{j} \right]_{\underline{x}}^{\bar{x}}. $$

Similarly, if radial basis functions $\psi_j(x) = e^{-\nu(x-c_j)^2}$ define the upper and lower bounds, the area metric is proportional to:

$$ \mathcal{A}(\theta) \propto \sum_{j=0}^{n_u} u_j \left[ \operatorname{erf}(\sqrt{\nu}(x - c_j)) \right]_{\underline{x}}^{\bar{x}} - \sum_{j=0}^{n_l} l_j \left[ \operatorname{erf}(\sqrt{\nu}(x - c_j)) \right]_{\underline{x}}^{\bar{x}}. $$

were erf(·) is the error function and $c_j$ are the knots values in $\mathcal{X}$ defining the centers of the radial basis $\psi_j(x)$.

Note that the area $\mathcal{A}(\theta)$ gives a measure of the accuracy (or precision) of an IPM as the average width of interval predictions. In this context, accuracy must not be confused with an IPM reliability, the probability of new samples falling within the predicted intervals. As example, consider an IPM defined by intervals of width 0.02 around a best fit function $\hat{f}(x)$. This IPM is considered very precise because it yields prediction intervals that are very narrow, i.e., $\hat{f}(x) \pm 0.01$. However, the reliability of the interval can be very small if the bounds do not quantify well the level of uncertainty affecting the process. For this reason, precision (or accuracy) and reliability (or probability of correct prediction) are competitive objectives and both are metrics of the goodness of an IPM.

## 2.4. IPM optimization programs

In this section we will review traditional IPM identification methods. The program in Eq. (3) presents an IPM identification method that minimizes the area metric $\mathcal{A}(\theta)$ while constraining all the available samples within $\mathcal{I}_N(x; \theta)$.

<table><caption>Table 1: A list of commonly adopted basis functions. The parameter $c_i$ are user-defined knots values in $\mathcal{X}$ and $\nu$ is a scale parameter defining the spread of a radial basis from its center $c_i$.</caption><thead><tr><th>Name</th><th>Basis function</th></tr></thead><tbody><tr><td>Polynomial basis</td><td>ψ<sub>j</sub>(x<sub>i</sub>) = x<sub>i</sub><sup>J</sup></td></tr><tr><td>Radial (Gaussian) basis</td><td>ψ<sub>j</sub>(x<sub>i</sub>) = exp(-ν||c<sub>j</sub> - x<sub>i</sub>||)</td></tr><tr><td>Fourier basis</td><td>ψ<sub>j</sub>(x<sub>i</sub>) = cos(c<sub>j</sub>x<sub>i</sub>) for i odd</td></tr><tr><td>Fourier basis</td><td>ψ<sub>j</sub>(x<sub>i</sub>) = sin(c<sub>j</sub>x<sub>i</sub>) for i even</td></tr><tr><td>Sigmoid/Soft Thresholds</td><td>ψ<sub>j</sub>(x<sub>i</sub>) = <sup>1</sup>&frasl;<sub>1+exp(-c<sub>j</sub>·x<sub>i</sub>)</sup></td></tr></tbody></table>
**IPM program 1. (Full data enclosure with no exceptions [17])**

$$ 
\begin{array}{ll}
\theta^* = \arg \min_{\theta} & \mathcal{A}(\theta) \\
\text{s.t.} & g(\delta_i; \theta) \le 0, \quad i = 1, \dots, N,
\end{array}
\qquad (3)
$$

where $\theta^*$ is the optimized vector of fitting coefficients (the IPM design), and the function $g(\delta; \theta)$ defines joint constraints on the upper and lower bounds.

The random samples in the data set $\mathcal{D}_N$ define $N$ constraints in Eq. (3). The constraints are given by:

$$ g(\delta_i; \theta) = \max \{g_u(\delta_i; \theta), g_l(\delta_i; \theta)\}, \quad i = 1, \dots, N, \qquad (4) $$

where

$$ g_u(\delta_i; \theta) = y_i - f_u(x_i; \theta), \quad i = 1, \dots, N, \qquad (5) $$

define the data-enclosure constraints for the upper bound $f_u$, and

$$ g_l(\delta_i; \theta) = f_l(x_i; \theta) - y_i, \quad i = 1, \dots, N, \qquad (6) $$

define the constraints for the lower bound. To ensure $f_u(x; \theta) > f_l(x; \theta)$ for all $x \in \mathcal{X}$, the following upper dominance constraint is included in the optimization program:

$$ g_d(\delta; \theta) = \sup_{x \in \mathcal{X}} [f_l(x; \theta) - f_u(x; \theta)] \le 0. \qquad (7) $$

If $\psi_j(x)$ is a polynomial basis, a sum-of-squares formulation of this constraint can be used, see e.g. [55] further details. In this work, the upper dominance constraint is imposed numerically by enforcing non-positivity of the maximum of the difference, $f_l(x; \theta) - f_u(x; \theta)$, evaluated on an arbitrary number of points $x_0, x_1, \dots, x_m \in \mathcal{X}$.

The IPM resulting from program (3) tightly encloses all the samples $(x_i, y_i) \in \mathcal{D}_N$. This is due to the $N$ constraints defined by (5) and (6), which must be satisfied with no exceptions for $\theta^*$ to be feasible. A drawback of (3) is the low accuracy (low precision due to the large width of the predictive bonds) when outliers are affecting the data set. In this case, the area metric can be further minimized by discarding some of the data points in $\mathcal{D}_N$. To this end, the authors of [2] proposed to optimally discard samples from $\mathcal{D}_N$ while minimizing $\mathcal{A}(\theta)$ and enclosing between the bounds the remaining samples [2]. The optimization program in Eq. (8) summarizes this approach.

**IPM program 2. (Optimal removal of $N_o$ samples [2])**

$$ 
\begin{array}{ll}
\langle \theta^* \rangle = \arg \min_{\theta} & \mathcal{A}(\theta) \\
\text{s.t.} & g(\delta_i; \theta) \le 0, \quad \forall \delta_i \in \mathcal{D}_{N-N_o},
\end{array}
\qquad (8)
$$

where $N_o$ is the number of discarded samples, the set

$$ \mathcal{D}_{N-N_o} = \{ \delta_i \in \mathcal{D}_N : \hat{F}_{g(\delta_i; \theta)} \le \frac{N - N_o}{N} \}, \qquad (9) $$

contains the remaining $N - N_o$ samples defining constraints on $f_u$ and $f_l$, and $\hat{F}_{g(\delta_i; \theta)}$ empirical Cumulative Distribution Function (CDF) of the constraint values $g$.

Equation (8) discards $N_o$ samples corresponding to the tails of the empirical CDF of the constraint function $g$, which is a metric defining the extent by which the data point contributes to the spread of the interval prediction.

An alternative method was recently introduced by [75] and named Minimax layer approach. This method, differently from (3) and (8), minimizes the maximum absolute distance between a regression function $f(x_i; \theta)$ and the samples $y_i$.

**IPM program 3. (Minimax layer [75])**

$$ \theta^* = \arg \min_{\theta} \max_{i \in \{1, \dots, N\}} |f(x_i; \theta) - y_i|, \qquad (10) $$

where $f(x_i; \theta) = \psi^T(x_i)\theta$ is a regression function defined by a linear combination of basis functions $\psi^T$ and fitting coefficients $\theta$ as presented in equations (1).
The Minimax layer IPM optimizes $f(x_i; \theta)$ by minimizing the $L_\infty$-norm be the errors [75], i.e., by minimizing the maximum of the distance between $f$ and the samples of the dependent quantity $y_i$, $i = 1, ..., N$. Once an optimal $f(x_i; \theta^*)$ is obtained, a constant width IPM is constructed as follows [75]:

$$ \mathcal{I}_N(x; \theta^*) = [f(x; \theta^*) - \gamma_i^*; f(x; \theta^*) + \gamma_i^*], $$

where $\gamma_i^* = |f(x_i; \theta^*) - y_i|$ is the $i^{th}$ highest empirical cost value, i.e., the $i^{th}$ largest difference between the regression function $f(x_i; \theta^*)$ and the samples of the dependent quantity. This IPM has a constant width of $2\gamma_i^*$ and the highest empirical cost value $\gamma_1^* = \max_{i \in \{1, ..., N\}} |f(x_i; \theta^*) - y_i|$ renders the IPM with the lowest accuracy (largest bounds) and highest reliability.

## 2.5. Convexity of the reviewed IPM programs

The functions $f_u$ and $f_l$ have a linear dependency with the optimization variables $\theta$. The area metric also displays a linear dependency in $\theta$. In fact, by assuming where $f_u(x; \theta) \ge f_l(x; \theta)$ for all $x \in \mathcal{X}$:

$$ \mathcal{A}(\theta) = \int_{\mathcal{X}} \left( \sum_{j=0}^{n_u} u_j \psi_j(x) - \sum_{j=0}^{n_l} l_j \psi_j(x) \right) dx = \sum_{j=0}^{n_u} u_j \int_{\mathcal{X}} \psi_j(x) dx - \sum_{j=0}^{n_l} l_j \int_{\mathcal{X}} \psi_j(x) dx. $$

The reviewed programs are convex with respect to the design parameters $\theta$ and can be solved very efficiently using a linear programming algorithm<fnref n="1" />. For further details on linear programming implementation of IPM see for instance [76].

## 2.6. Reliability of an IPM

The reliability of an IPM is the probability that $\mathcal{I}_N(x; \theta)$ will contain future samples from the DGM and is mathematically defined as follows:

$$ R(\theta) = \mathbb{P}[(x, y) \in \Delta : y \in \mathcal{I}_N(x; \theta)]. \quad (11) $$

Similarly, the probability of error (or failure probability) of an IPM is given by:

$$ P_f(\theta) = 1 - R(\theta) = \mathbb{P}[(x, y) \in \Delta : y \notin \mathcal{I}_N(x; \theta)]. \quad (12) $$

Note that the true values of $R(\theta)$ and $P_f(\theta)$ are unknown in practice. This is due to a limited availability of samples $(x, y) \in \mathcal{D}_N$ and, consequently, an epistemic uncertainty affecting the probability $\mathbb{P}$ generating the data. Sample-based approximations of $R(\theta)$ and $P_f(\theta)$ can be obtained from the available data as follows:

$$ \hat{R}(\theta) = \frac{\sum_{i=1}^{n} \mathbf{1}\{y_i \in \mathcal{I}_N(x_i; \theta)\}}{n}, \quad \hat{P}_f(\theta) = \frac{\sum_{i=1}^{n} \mathbf{1}\{y_i \notin \mathcal{I}_N(x_i; \theta)\}}{n}, \quad (13) $$

where $n$ is the number of new samples used for the estimation and $\mathbf{1}\{c_i\}$ is the indicator function for a condition $\{c_i\}$:

$$ \mathbf{1}\{c_i\} = \begin{cases} 1 & \text{if } c_i \text{ occurs} \\ 0 & \text{otherwise} \end{cases}. $$

Note that, by the weak law of large numbers, the errors between the estimators and the values of the true $R$ and $P_f$ converge asymptotically to an arbitrarily small value. As example for the reliability we have:

$$ \lim_{n \to \infty} \mathbb{P} \left[ |R(\theta) - \hat{R}(\theta)| \ge \epsilon \right] = 0, \quad (14) $$

where $\epsilon$ is an arbitrarily small margin of error. The number $n$ is always finite in practice and the difference between the true mean and the sample mean always non-zero.

The reliability is defined by a sequence Bernoulli trails where the binary outcome $\mathbf{1}\{c_i\} = 1$ has an (unknown) success chance $p$ and, thus, $R$ for $n$ samples follows a Binomial distribution, for which the following equality holds

$$ \mathbb{P}[R(\theta) - \hat{R}_m(\theta) \ge \epsilon] = \sum_{k=0}^{n(p-\epsilon)} \binom{n}{k} p^k (1-p)^{n-k}. $$

<footnote n="1">¹Note however, that the maximum operator preserves convexity but not linearity. To make the programs linear the constraints $g_u$ and $g_l$ should be enforced separately (imposing $2N$ constraints).</footnote>
The fact that $p$ is unknown does not allow the exact computation of this probability.

In the next sections, we will introduce Scenario optimization theory which allows quantifying the uncertainty on $R(\theta)$ and $P_f(\theta)$ by prescribing a formal mathematical expressions for the bound $\epsilon$. The bound quantifies the lack of data uncertainty which limits our knowledge of an IPM reliability. More in general, formal epistemic bounds can be obtained for the solution $\theta^*$ of a scenario optimization problem. Scenario optimization programs and formal reliability bounds prescribed via Scenario theory will be presented next.

### 3. A summary of Scenario theoretic epistemic reliability bounds

Consider a chance-constrained reliability program:

$$ \theta^* = \arg \min_{\theta \in \Theta} \{J(\theta) : R(\theta) > \alpha\}, \qquad (15) $$

where $\theta \in \Theta \subseteq \mathbb{R}^{n_\theta}$ is a vector of design variables constrained in a closed set $\Theta$, $J: \Theta \mapsto \mathbb{R}$ is a cost function and the reliability constraint $R(\theta) = \mathbb{P}[\delta \in \Delta : g(\delta; \theta) \le 0]$, enforces non-positiveness of $g$ with a probability at least $\alpha \in [0, 1]$. Problems like the ones defined by Eq. (15) are challenging to solve and often require approximating the probabilistic constraint, e.g., via Monte Carlo sampling, Most Probable Point method [77], Subset Simulation [78, 79], Line Sampling [80], Importance Sampling [81, 82] or other means [83, 84, 85]. The vast majority of the existing methods require a model of the uncertainty for $\delta$ to carry on with the analysis. However, this can lead to an overconfidence in the reliability $R(\theta^*)$ and to an under-estimation of the true level of uncertainty affecting the performance of $\theta^*$, especially under a lack of samples [6, 35].

Differently from traditional methods, Scenario programs gives a feasible solutions for Eq. (15) directly from $N$ samples of $\delta$ and without assumptions on the family of $\mathbb{P}[\cdot]$ generating the data. Scenario theory prescribes formal epistemic bounds on the reliability $R(\theta^*)$ which quantify the worst-case reliability among the possible $\mathbb{P}[\cdot]$ which are consistent with the observed samples. The bounds hold distribution-free, i.e., with no distributional assumptions on the underlying $\mathbb{P}[\cdot]$, and non-asymptotically, i.e., the inequality holds true for any number of sample $N$. In the next sections we will review scenario optimization programs, relevant definitions, assumptions, and theorems for bounding the reliability of scenario-based solutions.

#### 3.1. Scenario optimization programs

**Scenario program 1. (Hard-constrained, no exceptions)** A scenario optimization program with hard-constraints can be defined as follows:

$$ \theta^* = \arg \min_{\theta \in \Theta} \{J(\theta) : g(\delta_i; \theta) \le 0, i = 1, \dots, N, \} \qquad (16) $$

where $g(\delta; \theta) : \Theta \times \Delta \to \mathbb{R}$ is a constraint function defining $N$ samples constraints. The constraints in (16) are hard because an optimal design $\theta^*$ to be feasible must satisfy $g(\delta_i; \theta^*) \le 0$, for all $i = 1, \dots, N$, with no exception. Note that the IPM program in Eq. (3) is a particular type of hard-constrained scenario program where $g$ and $J$ are convex in $\theta$ and define, respectively, the bounds-enclosure constraints and the area metric objective function.

In many cases, (16) can lead to a feasibility issue if $\nexists \theta \in \Theta : g(\delta; \theta) \le 0$, for all $\delta \in \mathcal{D}_N$. To amend for this deficiency, a scenario program with discarded samples can be considered as presented by (2).

**Scenario program 2. (Hard-constrained, discarded samples)** A feasible sampling-and-discarding reformulation of program (16) is given by:

$$ \theta^* = \arg \min_{\theta \in \Theta} \{J(\theta) : g(\delta_i; \theta) \le 0, \delta_i \in \mathcal{D}_N \setminus \mathcal{O}\}, \qquad (17) $$

where $\mathcal{D}_N \setminus \mathcal{O}$ is the set difference between the data and a set of discarded samples $\mathcal{O}$. The cardinality of the set of discarded scenarios is $N_o = |\mathcal{O}|$.

The set $\mathcal{O}$ includes, for instance, samples which make (16) unfeasible, outliers, or scenarios which prevent the objective function from further improving. A solutions of (17), compared to a solution of program (16), enjoys an improved cost $J(\theta^*)$. However, this generally comes at the cost of weaker generalization guarantees. Note that IPM program (8) is a particular type of hard-constrained scenario program with discarded samples where $\mathcal{D}_N \setminus \mathcal{O}$ is defined by a condition on the tails of the constraint function $g$, see Eq. (9).

An alternative approach to enforce feasibility is the scenario approach with constraints relaxation recently investigated by [54].
**Scenario program 3. (Soft-constrained)** An feasible reformulation of (16) is given by:

$$ \langle \theta^*, \zeta^* \rangle = \arg \min_{\substack{\theta \in \Theta \\ \zeta \in \mathbb{R}_+^N}} \left\{ J(\theta) + \rho \sum_{i=1}^N \zeta_i : g(\delta_i; \theta) \le \zeta_i, \delta_i \in \mathcal{D}_N \right\}, \quad (18) $$

where $\zeta$ is a $N$-dimensional vector of non-negative slack variables and $\rho$ defines the cost of violations. A $\zeta_i = 0$ means that the hard-constraint imposed by the $i^{th}$ sample is satisfied, i.e., $g(\delta_i; \theta) \le 0$. On the other hand, a $\zeta_i > 0$ implies a violation of the hard-constraint. The scenario program defined by Eq. (18) always admits a feasible solution and if the cost parameter $\rho \to \infty$ the soft-constrained program goes back to the original formulation in Eq. (16).

A final scenario optimization program reviewed in this work is the following minimax optimization program. The epigraphic formulation of the problem is defined as follows:

**Scenario program 4. (minimax or $L_\infty$ minimization)**

$$ \langle \theta^*, \gamma^* \rangle = \arg \min_{\substack{\theta \in \Theta \\ \gamma \in \mathbb{R}}} \left\{ \gamma : \max_{i \in \{1, \dots, N\}} g(\delta_i; \theta) \le \gamma \right\}, \quad (19) $$

where $\gamma^*$ is the minimized $L_\infty$ norm (maximum) of $g(\delta_i; \theta^*)$, $i = 1, \dots, N$ at the optimum $\theta^*$.

The programs seek a design $\theta^*$ which minimizes the in-sample maximum of the constraint function $g$. The IPM Minimax layer program in Eq. (10) is special case of scenario program 4 where the constraint function is $g(\delta_i; \theta) = |f(x_i; \theta) - y_i|$ and is, therefore, convex in the design variables $\theta$. For earlier works on minimax robust approaches see, e.g., [86] and [49].

## 3.2. Definitions and assumptions for the derivation of scenario-based epistemic reliability bounds

**Definition 1. (Violation probability, or error probability)** The probability

$$ P_f(\theta^*) = \mathbb{P}[\delta \in \Delta : \theta^* \notin \Theta_\delta], \quad (20) $$

is called violation probability, or error probability. The quantity $\Theta_\delta$ is the feasibility set induced by $\delta$. Given a reliability parameter $\epsilon \in [0, 1]$ a solution $\theta^*$ is called $\epsilon$-robust if $P_f(\theta^*) < \epsilon$.

The quantity $P_f(\theta^*)$ gives the 'exact' probability that a solution $\theta^*$ will fail to meet the desired requirements, i.e., the probability that, given a new observation $\delta \in \Delta$, the constraint imposed by $g(\theta^*, \delta)$ in correspondence of $\theta^*$ will not be satisfied. The feasibility set induced by the constraints in program (1) is defined as follows:

$$ \Theta_\delta = \{\theta \in \Theta : g(\delta; \theta) \le 0\}. \quad (21) $$

Similarly, the feasibility set induced by the constraint in program (4) is given by:

$$ \Theta_\delta = \{\theta \in \Theta : g(\delta; \theta) \le \gamma^*\}. \quad (22) $$

Note that for IPM programs the violation probability corresponds to the probability of the interval failing to predict a future observation correctly.

**Definition 2. (Non-reducible support set)** A support set $S \subseteq \mathcal{D}_N$ is a k-tuple $S = \{\delta_{i_1}, \dots, \delta_{i_k}\}$ for which the solutions of a scenario program with constraints defined by $\mathcal{D}_N$ is identical to the same program with $S$ in place of $\mathcal{D}_N$. A set $S^* \subseteq S$ is non-reducible if, for any $\delta \in S^*$, the constraints defined by $S^* \setminus \delta$ leads to an optimum that differs from the one of the original problem. The support set $S^*$ is of minimal cardinality.

**Definition 3. (Complexity)** The dimension of $S^*$ denotes the complexity of a scenario solution and is given by $s_N^* = |S^*|$, where $| \cdot |$ is the cardinality operator. Note that scenario program generally admits several support sets and a the set of minimal cardinality $S^*$ renders the tightest bounds.

**Assumption 1. (Existence and uniqueness)** For any positive integer $N \in \mathbb{N}_0$ and scenario set $\mathcal{D}_N$, the optimal solution $\theta^*$ of a scenario program exists and is unique.

Existence can be enforced by constraining the optimization on a compact domain $\Theta$, this avoids problems which may arise when the objective function decreases when $\theta$ drifts to infinity. Concerning the uniqueness of the solution, if multiple optimal solutions exist in $\Theta$ a tie-break rule can be implemented, e.g. selecting the solution with minimum norm or optimizing additional convex functions in $\theta$.
**Assumption 2.** (Non-degeneracy) For any positive integer $N \in \mathbb{N}_0$ and scenario set $\mathcal{D}_N$, the solution of the scenario program with $N$ constraints coincides with probability 1 with the solution of the same program with constraints defined by $S$.

When (16) is convex, non-degeneracy is a mild assumption since support constraints are always active constraints (but the converse not always remains true). In the general non-convex case however, $S$ might include non-active constraints. In this case, the removal of a non-active scenario constraint, e.g., $g(\delta_i, \theta^*) < 0$, yield a new optimal solution having a smaller cost $J(\theta^*)$ [36].

**Assumption 3.** (Non-accumulation) For every $\theta \in \Theta$ the statement $\mathbb{P}[\delta : g(\delta; \theta) = a] = 0$ for any $a \in \mathbb{R}$ holds true.

This is equivalent to assume that the function $g(\delta; \theta)$ does not have concentrated mass. This assumption is generally satisfied when $\delta$ admits a probability density function.

### 3.3. An overview of epistemic reliability bounds

Scenario theory gives formal epistemic bounds on the probability of an optimized scenario design $\theta^*$ failing to comply with new observations, i.e., formal bounds on $P_f(\theta^*)$ and $R(\theta^*)$. These bounds are also known as “generalization” bounds, in the sense that a bound gives guarantees on the solution’s ability to perform well when future not yet observed data is collected. Note that the complexity of a scenario decision is a fundamental concept in the derivation of the generalization error bounds. In fact, for an equal number of available scenarios $N$, solutions having high complexity $s_N^*$ yield weaker reliability guarantees. Also, scenario-based bounds are ‘epistemic bounds’ because the uncertainty on the reliability score can be reduced by collecting more samples and scenario-bounds gets tighter when the data set size $N$ increases. In the next sections we will review Scenario theoretic bounds for the reviewed scenario programs.

### 3.4. Bounds for hard-constrained convex problems

**Theorem 1.** [24, Theorems 2.4] Consider a convex scenario program defined as in (16). Under assumptions 1, 2, stationary $\mathbb{P}$ and IID samples in $\mathcal{D}_N$ and $N > n_\theta$ the distribution of $P_f(\theta^*)$ is bounded by a Beta distribution:

$$ \mathbb{P}^N[P_f(\theta^*) \le \epsilon] \ge 1 - \sum_{k=0}^{n_\theta-1} \binom{N}{k} \epsilon^k (1-\epsilon)^{N-k} = 1-\beta, \quad (23) $$

where $n_\theta$ is the number of design variables and $\beta \in [0, 1]$ is a confidence parameter.

The parameter $\epsilon$ gives an upper bound on the probability of error and selecting a confidence parameter $\beta$, it can be expressed as a function $\bar{\epsilon}(N, \beta, n_\theta)$ of the number of design variables $n_\theta$, the number of samples $N$, and a confidence level $\beta$ by solving the right-hand side of Eq.(23). This bound quantifies the epistemic uncertainty affecting $P_f$ and is obtained before observing the solution $\theta^*$ of Eq. (16), i.e., it is a-priori computed. Theorem 1 was derived from the Helly’s theorem and showing that $s_N^* \le n_\theta$ for hard-constrained convex scenario program. In other words, the complexity of the solution of a convex scenario program is always capped by $n_\theta$ and the solution of the problem always admits a compressed representation of dimension $n_\theta$. For optimization program like in Eq. (16), when fully-supported, we have that $s_N^* = n_\theta$ with probability 1. For this type of problems, Eq. (23) has been proven to hold tightly, i.e., with the equality sign.

However, the majority of convex problems are only partially-supported and the solution $\theta^*$ can be often reconstructed by a number of constraints $s_N^* < n_\theta$. For these problems, Eq. (23) gives a conservative a-priori bound. A tighter bound can be obtained by an a-posteriori evaluation of the solution’s complexity. Theorem 2 presents the a-posteriori bound for convex hard-constrained scenario programs, also known as the wait-and-judge method [25].

**Theorem 2.** [25, Theorem 2] Consider a convex scenario program defined as in (16). Under assumptions 1, 2, stationary $\mathbb{P}$ and IID samples in $\mathcal{D}_N$ the solution $\theta^*$ of (16) satisfies

$$ \mathbb{P}^N[P_f(\theta^*) \le \epsilon(s_N^*)] \ge 1 - \beta, \quad (24) $$

where the reliability $\epsilon(k) = 1 - t(k)$ is the unique solution in $[0,1]$ of the following polynomial equation in $t$ for any $k = 0, ..., n_\theta$:

$$ \frac{\beta}{N+1} \sum_{j=k}^{N} \mathfrak{B}_j(t; k) - \binom{N}{k} t^{N-k} = 0. \quad (25) $$

Here $\mathfrak{B}_j(t; k) = \binom{j}{k} t^{j-k}$ is a binomial expansion.
Theorem 2 gives a generalization bound $P_f(\theta^*) \le \epsilon(s_N^*)$ which holds up to a confidence level $1-\beta$. The value of the reliability parameter $\epsilon(s_N^*)$ is determined from (25) after observing $s_N^* = k$ support constraints which is always $\le n_\theta$ for convex problems. In contrast to the bound in (23), $\epsilon(s_N^*)$ is a-posteriori computed after enumeration of the set of support scenarios and its cardinality $s_N^*$ in correspondence of $\theta^*$. An important consequence of this theorem is that a-posteriori observing a solution that has complexity of $s_N^*$ form a dimension $n_\theta > s_N^*$ does not give equal robustness guarantees as obtaining a solution from a problem in dimension $n_\theta = s_N$ whose solutions are always representable by $s_N$ observations. That is, simple questions are more guaranteed than simple answers. The interested reader is reminded to [25] for a detailed discussion about this topic. Theorem 2 was proven by a characterization of Eq. (16) in terms of a generalized moment problem and duality theory based on this characterization.

### 3.5. Bounds for convex problems with discarded samples

An extension of Theorem 1 allows for $k$ samples $\delta$ to be intentionally discarded from $D_N$, for instance, the ones making program (16) unfeasible. As a result, the optimized $\theta^*$ enjoys an improved $J(\theta^*)$ for the cost of a weaker certificate of generalization [24].

**Theorem 3.** [87, Theorem 2.1] Consider a convex scenario program defined as in Eq. (17) where $k = N_0$ scenarios have been intentionally removed from the data set $D_N$. Under assumptions 1, 2 and stationary $\mathbb{P}$ and IID samples in $D_{N-k}$ the distribution of $P_f(\theta^*)$ is bounded by:

$$ \mathbb{P}^m[P_f(\theta^*) \le \epsilon] \ge 1 - \binom{k+n_\theta-1}{k} \sum_{i=0}^{k+n_\theta-1} \binom{m}{i} \epsilon^i (1-\epsilon)^{m-i} = 1-\beta, \quad (26) $$

where $n_\theta$ is the number of optimization variables and $\beta \in [0, 1]$ is a confidence parameter.

Theorem 3 is used in practice as follows: a desired confidence parameter $\beta$ and $N_0$ samples which defines the discarded set of scenarios $O$ are selected. Given the number of available scenarios $N$ and dimension of the design variables $n_\theta$ a reliability parameter $\epsilon(N, \beta, n_\theta, N_0)$ is computed from the right hand side of (26). This gives an upper bound $P_f(\theta^*) \le \bar{\epsilon}(N, \beta, n_\theta, N_0)$ on the violation probability which is a-priori computed, i.e., is obtained before observing $\theta^*$ solution of program 2.

### 3.6. Bounds for soft-constrained convex problems

To extend the scope of scenario-based generalization bounds to the soft-constrained case (18), the technical assumption of non-accumulation is required. For these problems the following theorem holds.

**Theorem 4.** [54, Theorem 4] Consider a convex scenario program as in Eq. (18). Given assumptions 1, 2, 3, stationary $\mathbb{P}$ and IID samples in $D_N$ the probability $P_f(\theta^*)$ is bounded by:

$$ \mathbb{P}^N[\underline{\epsilon}(k) \le P_f(\theta^*) \le \bar{\epsilon}(k)] \ge 1 - \beta, \quad (27) $$

where $k = s_N^*$ is the number of support constraints and $\underline{\epsilon}(k) = \max\{0, 1 - \underline{t}(k)\}$, $\bar{\epsilon}(k) = 1 - \bar{t}(k)$ and $\{\underline{t}(k), \bar{t}(k)\}$ are solutions of a polynomial equation in $t$:

$$ \mathfrak{B}_N(t; k) = \frac{\beta}{2N} \sum_{j=k}^{N-1} \mathfrak{B}_j(t; k) + \frac{\beta}{6N} \sum_{j=N+1}^{4N} \mathfrak{B}_j(t; k), \quad (28) $$

where $k = 1, \dots, N-1$ is the number of support constraints. For the special case $k = N$, the upper bound is set to $\bar{\epsilon}(k) = 1$ and the lower bound is obtained solving

$$ 1 = \frac{\beta}{6N} \sum_{j=N+1}^{4N} \mathfrak{B}_j(t; k). \quad (29) $$

Theorem 4 gives a-posteriori upper and lower bounds on the violation probability of the solution of (18). The bounds can be computed in practice from the number of constraints $N$ by selecting a confidence $\beta$ and enumerating the support constraints $s_N^*$ i.e., $[\underline{\epsilon}(N, \beta, s_N^*), \bar{\epsilon}(N, \beta, s_N^*)]$. Note that, for soft-constrained convex programs, $s_N^* = \sum_{i=1}^N 1(g(\delta_i; \theta^*) \ge 0)$. The number of support constraints coincides with the number of scenarios for which $g(\delta_i; \theta^*) \ge 0$ at the optimum $\theta^*$.
### 3.7. Bounds for minimax problems

Let us define the ordered statistic of the constraint function values as $\gamma_1^* \ge \gamma_2^* \ge \dots \ge \gamma_N^*$ where the quantity $\gamma_i^* = g(\delta_i; \theta^*)$ is the $i^{th}$ higher value of the constraint $g$ in correspondence of the optimum $\theta^*$. The probability of error for an empirical level $\gamma_i^*$ is given by $P_{f,i} = \mathbb{P}[\delta \in \Delta : g(\delta; \theta^*) > \gamma_i^*]$, and the corresponding reliability is $R_i = \mathbb{P}[\delta \in \Delta : g(\delta; \theta^*) \le \gamma_i^*]$. Note that $\gamma_i \ge \gamma_j$ for any $i < j$ and, thus, the empirical probabilities of error are $P_{f,i} \le P_{f,j}$ for $i < j$. Theorem 5 prescribes a formal bound on the joint distribution of the probabilities of error $P_{f,n_{\theta}+1}, \dots, P_{f,N}$.

**Theorem 5.** [31, Theorems 1-2] For convex minimax program as in Eq. (4), the joint probability distribution function of $P_{f,n_{\theta}+1}, \dots, P_{f,N}$ is bounded by

$$ \mathbb{P}^N[P_{f,i} \le \epsilon_{n_{\theta}+1}, \dots, P_{f,N} \le \epsilon_N] \le F_{\theta,N}(\epsilon_{n_{\theta}+1}, \dots, \epsilon_N) = \beta, \quad (30) $$

where

$$ F_{n_{\theta}, N} = \frac{N!}{n_{\theta}!} \int_{0}^{\epsilon_{n_{\theta}+1}} \tau^{n_{\theta}}_{n_{\theta}+1} \dots \int_{0}^{\epsilon_N} \mathbf{I}_{\{0 \le \tau_{n_{\theta}+1} \le \dots \tau_N \le 1\}} d\tau_N \dots d\tau_{n_{\theta}+1}, $$

is the ordered Dirichlet distribution bounding the joint probability distribution function of the empirical errors. For non degenerate problems, the bound (30) is tight and holds with equality sign.

Since $P_{f,n_{\theta}+1} \ge P_{f,i}$ for $i = 1, \dots, n_{\theta}$, the bound $\epsilon_{n_{\theta}+1}$ on the level $P_{f,n_{\theta}+1}$ automatically applies also to $P_{f,i}$, $i = 1, \dots, n_{\theta}$.

In practice, Theorem 5 can be used to construct an upper and lower bound on the empirical CDF of the probability of a solution violating an empirical level $\gamma$. This is a non-parametric probability box which bounds the reliability at the empirical violation levels $\gamma$. For an overview on probability boxes the interested reviewer is reminded to [9, 88, 89, 90, 91] and propagation strategies on different applications. The upper bound of the probability box is constructed as follows: [31]:

$$ \bar{F}(\gamma) = \begin{cases} 1 & \text{if } \gamma \ge \gamma_{n_{\theta}+1} \\ 1 - \epsilon_i & \text{if } \gamma_{i+1} < \gamma \le \gamma_i, (i = n_{\theta}+1, \dots, N-1), \\ 1 - \epsilon_N & \text{if } \gamma < \gamma_N \end{cases}, \quad (31) $$

whilst the lower bound is given by,

$$ \underline{F}(\gamma) = \begin{cases} 1 - \bar{\epsilon}_1 & \text{if } \gamma \ge \gamma_1, \\ 1 - \bar{\epsilon}_i & \text{if } \gamma_i \le \gamma \le \gamma_{i-1}, (i = 2, \dots, N), \\ 0 & \text{if } \gamma < \gamma_N \end{cases}, \quad (32) $$

where $\bar{\epsilon}_i$ and $\epsilon_i$ are obtained for a confidence level $\beta$ from Theorem 4. The function $\bar{F}(\gamma)$ and $\underline{F}(\gamma)$ are the upper and lower CDFs bounding empirical reliability levels

$$ \underline{F}(\gamma) \le \mathbb{P}[\delta \in \Delta : g(\delta; \theta^*) \le \gamma] \le \bar{F}(\gamma), $$

with a confidence at least $1 - 2\beta$.

### 3.8. Bounds for non convex problems

**Theorem 6.** [28, Theorems 1, Equation (7)] Consider a non-convex scenario program defined like in Eq. (16). Under assumptions 1, 2, and IID samples in $D_N$, the distribution of the violation probability $P_f(\theta^*)$ is bounded by,

$$ \epsilon(N, \beta, s_N^*) = \begin{cases} 1 & \text{if } s_N^* = N \\ 1 - \left(\frac{\beta}{N \binom{N}{N}}\right)^{\frac{1}{N-1}} & \text{otherwise,} \end{cases} \quad (33) $$

where $s_N^*$ is the number of support scenarios and $\beta \in [0, 1]$ is the selected confidence.

Note that this theorem offers a conservative upper bounds on the violation probability which holds for non convex programs defined as in (16), i.e., programs for which $g(\delta; \theta)$ or $J(\theta)$ are non-convex functions of $\theta$. If an optimization program defined as in Eq. (16) is convex, the use of less conservative a-posteriori Wait-and-judge bounds as defined in theorem (2) is advisable. For the remainder of this paper we will focus on convex programs and reliability bounds for convex problems.
### 3.9. Summary of the reviewed bounds

Scenario theory allows deriving powerful reliability bounds for the solutions of scenario optimization programs. The number of samples $N$, a measure of the complexity of a solution, and a confidence parameter $\beta$ determine the bound value. In particular, the complexity of a solution can be defined as the number of scenarios that are needed to reproduce the optimum. For convex problems, and given some technical assumptions, the complexity of a solution is always upper bounded by the number of optimization variables. The unique solution of convex optimization with $N$ scenario constraints can be always reconstructed by solving the same optimization problem with maximum $n_{\theta}$ samples constraints, e.g., the ones that are active at the optimum. Hence, for convex problems scenario-based reliability bound can be computed from the number of optimization variables $n_{\theta}$ before solving the optimization and observing the optimized solution, i.e., a-priori [24, 87]. However, a-priori bounds may result in over-conservative wide bounds when the complexity $s_N^*$ is much lower than the number of decision variable $n_{\theta}$.

In contrast, a-posteriori evaluation of the solution's complexity can provide a tighter quantification of the episodic uncertainty affecting its reliability [25, 92, 54]. A-posteriori bounds requires computing the number support constraints $s_N^*$ that is measure of the complexity of the solution. Note that, in contrast to $n_{\theta}$, the number $s_N^*$ is only available after optimizing $\theta^*$ and is a random variable that depends on the multi-samples extraction $\mathcal{D}_N$. Table 2 summarizes the reviewed bounds and theorems for convex and non-convex problems and computed a-priori or a-posteriori. The check mark (X) indicates that, to the best author knowledge, scenario bounds are not yet available or not explicitly defined. For instance, this is the case for a-priori bounds for non-convex scenario programs and for convex soft-constrained problems. The unavailability of a-priori bounds for these types of problems is due to our inability to a-priori bound the complexity of their solutions, i.e., in principle all the available samples can be indispensable to replicate the solution and are therefore supports. In contrast, the complexity of the solutions of hard-constrained problems with discarded samples can be obtained a-priori and a-posteriori from the dimension $n_{\theta}$ (or the number of support scenarios $s_N^*$) and the number of discarded samples $N_o$. Figure 1 displays the bounds (in log-scale) for increasing $N$. The bounds are obtained for a confidence parameter $\beta = 10^{-6}$, assuming $n_{\theta} = 100$, $s_N^* = 95$ and $N_o = 5$ discarded samples. For additional details and a numerical comparison of the bounds see, for instance, [93].

<table><caption>Table 2: Summary of the reviewed reliability bounds $\epsilon$ for different scenario programs, dependency with relevant factors in the problems, reference theorem and the bounds a-priori or a-posteriori computability. The check mark (X) means that, to the best author knowledge, a-priori or a-posteriori scenario bounds on the violation probability are not available or not explicitly defined.</caption><thead><tr><th rowspan="2">Type</th><th colspan="2">Scenario</th><th colspan="2">convex</th><th colspan="2">non-convex</th><th colspan="2">Type</th><th rowspan="2">Ref.</th></tr><tr><th>Program</th><th>Bounds</th><th>Theorem</th><th>Bounds</th><th>Theorem</th><th>a-prio.</th><th>a-post.</th></tr></thead><tbody><tr><td>Hard-constrained</td><td>1</td><td>ε̄(N, β, n<sub>θ</sub>)</td><td>1</td><td>X</td><td>X</td><td>✓</td><td></td><td>[24]</td></tr><tr><td>Hard-constrained</td><td>1</td><td>ε̄(N, β, s<sub>N</sub><sup>*</sup>)</td><td>2</td><td>ε̄(N, β, s<sub>N</sub><sup>*</sup>)</td><td>6</td><td></td><td>✓</td><td>[25] [28]</td></tr><tr><td>Hard-constrained (discarded)</td><td>2</td><td>ε̄(N, β, n<sub>θ</sub>, N<sub>o</sub>)</td><td>3</td><td>X</td><td>X</td><td>✓</td><td></td><td>[87]</td></tr><tr><td>Hard-constrained (discarded)</td><td>2</td><td>ε̄(N, β, s<sub>N</sub><sup>*</sup> + N<sub>o</sub>)</td><td>2</td><td>X</td><td>X</td><td>✓</td><td></td><td>[25]</td></tr><tr><td>Soft-constrained</td><td>3</td><td>X</td><td>X</td><td>X</td><td>X</td><td>✓</td><td></td><td>[]</td></tr><tr><td>Soft-constrained</td><td>3</td><td>ε̄(N, β, s<sub>N</sub><sup>*</sup>), ε̄(N, β, s<sub>N</sub><sup>*</sup>)</td><td>4</td><td>X</td><td>X</td><td>✓</td><td></td><td>[54]</td></tr><tr><td>Minimax</td><td>4</td><td>ε̄(N, β, n<sub>θ</sub>)</td><td>1</td><td>X</td><td>X</td><td>✓</td><td></td><td>[24] [94]</td></tr><tr><td>Minimax</td><td>4</td><td>ε<sub>1</sub>(N, β, n<sub>θ</sub>), ..., ε<sub>N</sub>(N, β, n<sub>θ</sub>)</td><td>5</td><td>ε̄(N, β, s<sub>N</sub><sup>*</sup>)</td><td>6</td><td></td><td>✓</td><td>[31] [28]</td></tr></tbody></table>

Table 2 presents the bounds as functions $\epsilon = \epsilon(N, \beta, ...)$, which describes the most typical problems of bounding the reliability of a scenario solution given the number of scenarios a confidence, the number of discarded samples, etc. As example, Finally, note that the theorems presented in this work can be numerically inverted to help in the solution of other types of problems, for instance [17]:

* *Design of an identification experiment, N = N(ε, β, n<sub>θ</sub>, N<sub>o</sub>).* The problem of finding the minimum number of observations N which are necessary to obtain a certain reliability score $\epsilon$ given $\beta$, $n_{\theta}$ and $N_o$.

* *Maximal complexity and model selection, n<sub>θ</sub> = n<sub>θ</sub>(N, ε, β, N<sub>o</sub>).* The problem of determining the maximum number of the fitting coefficients $n_{\theta}$ to achieve a reliability at least $\epsilon$ given N samples, $N_o$ outliers, and a confidence $\beta$. As example, consider a problem where we target a reliability of at least 0.99 at a confidence $1 - 10^{-6}$ and $N = 5000$ scenarios are available to optimize a decision (of which $N_o = 5$ are removed). If we define our optimization problem to only include $n_{\theta} \le 7$ design variables, we can ensure a-priori a reliability level of at least 0.99.

* *Data cleansing, N<sub>o</sub> = N<sub>o</sub>(N, ε, β, n<sub>θ</sub>).* The process of detecting and removing corrupt or inaccurate data. Can help answer questions such as “How many data points can we discard while retaining a reliability of at least $\epsilon$?”
<figure><img src="image_1.png" /><figcaption>Figure 1: A comparison of the scenario bounds summarized in Table 2.</figcaption></figure>

## 4. The proposed soft-contained IPM with certified lower and upper epistemic reliability bounds

In contrast to traditional hard-constrained approaches, a soft-constrained IPM program softens the hard-constraints $g(\delta; \theta) \le 0$ on the lower and upper bounds. A relaxation of the constraints enables the condition $g(\delta; \theta) \le 0$ to be violated, i.e., some of the data points can fall outside from the optimized bounds if this yields a significant improvement of the area metric. Consider the following soft-constrained IPM training program.

**IPM program 4.** (Soft-Constrained IPM design)

$$ \begin{array}{ll} \langle \theta^*, \zeta^* \rangle = \arg \min_{\substack{\theta, \\ \zeta_i \ge 0}} & \mathcal{A}(\theta) + \rho \sum_{i=1}^{N} \zeta_i : \\ \text{s.t.} & f_i(x_i; \theta) - y_i \le \zeta_i, \quad i = 1, \dots, N, \\ & y_i - f_u(x_i; \theta) \le \zeta_i, \quad i = 1, \dots, N, \end{array} \qquad (34) $$

where $\zeta^*$ is a $N$-dimensional vector of non-negative slack variables relaxing the samples constraints on the upper and lower bounds.

If $\zeta_i^* > 0$ at the optimum, the sample $(x_i, y_i)$ will fall outside the bounds defined by $\mathcal{I}_N(x; \theta^*)$. Conversely, if $\zeta_i^* = 0$ the original hard-constraint on the lower and upper bounds are satisfied and the sample $(x_i, y_i)$ falls within the predicted interval, i.e., $y_i \in \mathcal{I}_N(x_i; \theta^*)$. Note that if the cost parameter $\rho \to \infty$ the soft-constrained program goes back to the original formulation (3) where all the sample must fall within the predictive bounds. Softening of hard-constraints yields IPMs designs with a lower value of the optimized area $\mathcal{A}(\theta^*)$, i.e., the resulting accuracy will be better than the hard-constrained counterpart. This is due to the (in average) narrower predictive interval. However, this generally leads to weaker guarantees on the reliability of the predictor.

As discussed, the complexity of convex programs can be generally bounded a-priori. However, the complexity of the proposed soft-constrained method cannot be a-priori bounded. This is due to the additional slack variables $\zeta$ leading to a decision space of dimension $N + n_\theta > N$ which can therefore lead to a solution having complexity $s_N^* > n_\theta$. Nonetheless, the reliability of a solution of (34) can be bounded a-posteriori using Theorem 4 which give both upper and lower bounds on the error probability of a soft-constrained IPM as follows:

$$ \mathbb{P}^N [\epsilon \le \mathbb{P}[(x, y) \in \Delta : y \notin \mathcal{I}_N(x; \theta^*)] \le \bar{\epsilon}] \ge 1 - \beta. $$

By selecting a confidence level, an epistemic interval on the probability of error is derived as follows:

$$ P_f(\theta^*) \in [\underline{\epsilon}(N, \beta, s_N^*), \bar{\epsilon}(N, \beta, s_N^*)]. $$
Similarly, an interval on the IPM reliability holds with the same confidence level and is given by:

$$R(\theta^*) \in [1 - \bar{\epsilon}(N, \beta, s_N^*), 1 - \underline{\epsilon}(N, \beta, s_N^*)].$$

The program in Eq. (34) is convex in $\theta$ and, thus, $s_N^*$ can be easily obtained from the number of violating scenarios (for which $\zeta_i^* > 0$) and the number of support scenarios on the bounds, scenario that if removed from the set of constraints improve the value of the objective function $\mathcal{A}$. For convex programs, a support constraint is always active at the optimum<fnref n="2" /> and a conservative estimate of $s_N^*$ is given by $s_N^* = |\mathcal{S}|$ where the set of support constraints is

$$\mathcal{S} = \{(x_i, y_i) \in \mathcal{D}_N : \zeta_i^* > 0, \lor y_i - f_u(x_i; \theta^*) = 0, \lor f_l(x_i; \theta^*) - y_i = 0\},$$

where $\lor$ is the logic disjunction defining the ‘or’ boolean operator.

## 5. Case studies

In the following section we will compare the accuracy of IPMs identified via the optimization methods described in Sections 2 and 4. The reliability bounds presented in Section 3.3 will be used to assess the generalization of the predictors and their validity will be verified empirically. The following examples will be considered:

Case study I: a trigonometric function affected by a noise term.

Case study II: modeling the response of dynamic system affected by uncertainty.

Case study III: a frequency response function of a car suspension arm affected by a structural damage with uncertain property, i.e., a cracks of unknown length.

The DGMs are controlled by the analyst and can be used to estimate the ‘true’ IPM reliability and assess the validity and tightness of the generalization error. For fairness of comparison, we solve all the IPM programs using the MATLAB ‘fmincon’ solver equipped with the sequential quadratic programming ‘sqp’ algorithm. The solver runs on an Intel(R) Core(TM) i7-8750H CPU machine, with a 2.2 GHz processor 32 Gb RAM.

### 5.1. Case study I

The DGM of the first case study is defined as follows [2, 55]:

$$y(x; \sigma) = x^2 \cos(x) - e^{-x^2} \sin(3x) - x - \cos(x^2) + x\mathcal{N}(0, \sigma) \quad (35)$$

where $x \in \mathcal{X}$ is an input variable in $\mathcal{X} = [-10, 5]$, and $\mathcal{N}(0, \sigma)$ is a Gaussian distributed random noise with zero mean and standard deviation $\sigma$. A data set $\mathcal{D}_N$ is obtained by uniform sampling of $N$ realization of $x_i$ in $\mathcal{X}$ and computing $y_i$ with a randomized noise value.

#### 5.1.1. Results for hard-constrained IPM programs: a-posteriori vs a-priori bounds

A data set $\mathcal{D}_{N=150}$ is sampled from Eq. (35) with $\sigma = 1$ and an optimal IPM design is determined solving program (3). We define the functions $[f_l, f_u]$ as linear combinations of polynomial basis $\psi(x)$ and $n_l = n_u$ fitting coefficients ranging between 3 and 20. Thus, the number of optimization variables $n_\theta$ is between 6 and 40. Figure 2 displays two IPMs optimized selecting polynomials of order $n_l = 5$ (black dashed lines) and order $n_l = 11$ (red solid lines). Note that a total of sixteen support scenarios are preventing the optimization program from further reducing the area metric of the IPM of order $n_l = 11$ and are presented by green round markers in the Figure. Table 3 presents the numerical results. As figures of merit, we compare the computational cost of the optimization, the (normalized) IPM area, and the a-priori and a-posteriori reliability bounds. Theorems 1 and 2 compute the a-priori and a-posteriori bounds on the error probability. For the a-posteriori bounds, the complexity is the dimension of the design space $n_\theta$ whilst for the a-posteriori bounds the complexity of the IPM solution is determined by the random number of support constraints $s_{150}^*$ which is at most equal to $n_\theta$. Monte Carlo estimators of the true probability of error $\hat{P}_f$ are obtained from $10^5$ samples always lay within the bounds and presented in the last row of Table 3. A higher value of $n_\theta$ leads to more accurate predictions (lower area) because $[f_l, f_u]$ can more tightly fit the data if more coefficients are used. However, this comes for the cost of weaker generalization error bound. Consider for instance the IPMs $n_l = 5$ and $n_\theta = 11$ presented in Figure 2, which render $\mathcal{A} = 0.79$ and $\mathcal{A} = 0.42$, respectively. Although $n_l = 5$ is less accurate its is also guaranteed to be more reliable. In fact, the a-posteriori error bounds result

<footnote n="2">²Note that an active constraint is not always a support constraint and therefore $s_N^*$ computed as the number of active constraints can be larger than the support set of minimal cardinality. This over-estimation will result in a conservative estimate of the reliability bound.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: IPM bounds designed with $n_l = n_u = 5$ (dashed black lines) and $n_l = n_u = 11$ (red solid lines). The 150 samples are displayed by blue markers whilst green markers shows the support constraints laying on the bounds of degree 11.</figcaption></figure>

<table><caption>Table 3: Area metrics, computational time and number of support constraints for $N = 150$ and IPM trained with program (3). The a-posteriori and a-priori Scenario bounds from Theorems 1 and 2, are also presented for a confidence parameter $\beta = 10^{-6}$ (almost certainty) and verified via Monte Carlo (MC) estimation of the true failure probability.</caption><tbody><tr><td>n<sub>l</sub> = n<sub>u</sub></td><td>3</td><td>4</td><td>5</td><td>7</td><td>9</td><td>11</td><td>15</td><td>20</td></tr><tr><td>CPU Time [s]</td><td>0.02</td><td>0.1</td><td>0.3</td><td>1.7</td><td>6.3</td><td>17.1</td><td>49.5</td><td>111.7</td></tr><tr><td colspan="9">Accuracy</td></tr><tr><td>A(θ<sup>*</sup>)</td><td>0.85</td><td>0.84</td><td>0.79</td><td>0.66</td><td>0.54</td><td>0.42</td><td>0.33</td><td>0.32</td></tr><tr><td colspan="9">Epistemic bounds on P<sub>f</sub></td></tr><tr><td>n<sub>θ</sub></td><td>6</td><td>8</td><td>10</td><td>14</td><td>18</td><td>22</td><td>30</td><td>40</td></tr><tr><td>ε̄(n<sub>θ</sub>, N, β) a-priori</td><td>0.158</td><td>0.18</td><td>0.20</td><td>0.24</td><td>0.28</td><td>0.31</td><td>0.37</td><td>0.45</td></tr><tr><td>s<sub>N</sub><sup>*</sup></td><td>6</td><td>8</td><td>10</td><td>14</td><td>15</td><td>16</td><td>19</td><td>19</td></tr><tr><td>ε̄(s<sub>N</sub><sup>*</sup>, N, β) a-posteriori</td><td>0.19</td><td>0.21</td><td>0.23</td><td>0.27</td><td>0.28</td><td>0.29</td><td>0.32</td><td>0.32</td></tr><tr><td colspan="9">Validation</td></tr><tr><td>P̂<sub>f</sub>(θ<sup>*</sup>), MC</td><td>0.06</td><td>0.072</td><td>0.065</td><td>0.105</td><td>0.151</td><td>0.153</td><td>0.215</td><td>0.230</td></tr></tbody></table>

$\bar{\epsilon}(10, 150, \beta) = 0.23$ for $n_l = 5$ and $\bar{\epsilon}(16, 150, \beta) = 0.29$ for $n_l = 11$. In other words, the first IPM is guaranteed to predict correctly at least 77% of the future data whilst the IPM with $n_l = 11$, although more accurate, is only guaranteed to be 71% reliable. Note however that for large $n_\theta$ the a-posteriori complexity is generally $s_{150}^* < n_\theta$, and sometimes much lower than $n_\theta$. For these cases, an a-posteriori evaluation of the IPM reliability renders the stronger guarantees. We repeated the analysis with $N = 2000$ scenarios and $n_l = n_u = 20$. The first column of Table 4 presents the normalized accuracy and reliability results and a comparison between a-priori and a-posteriori bounds.

### 5.1.2. Outliers removal

The program defined by Eq. (8) is tested to optimally remove outlier when constructing the IPM. Polynomial basis of degree $n_l = n_u = 9$ are considered and 10% of the samples in $D_{150}$ are discarded ($N_o = 15$). Figure 3 compares an IPM obtained without discarding samples (dashed black lines) with an IPM resulting from program (8) with $N_o = 15$ discarded scenarios (round marker). The latter IPM led to a reduction in the area metric of about 25% compared to the case with $N_o = 0$, thus offering (in average) narrower predictive interval and higher predictions. However, this comes at the price of a weaker generalization error guarantee. The IPM without discarded scenarios results in a $1 - \epsilon = 0.72$ (see Table 3) whilst the IPM obtained discarding samples is only $R(\theta^*) > 1 - \epsilon = 0.46$ reliable. Validation samples successfully predicted by the IPM are presented in the bottom panel in Figure 3 by blue markers and violating scenarios by red markers. Table 4 presents the results for IPM program in Eq. (8) and for a larger data set containing $N = 2000$ scenarios and increasing the number of discarded samples.
<figure><img src="image_3.png" /><figcaption>Figure 3: In the top panel: bounds resulting program (8) with and $N_o = 0$ (dashed lines) and $N_o = 15$ (solid lines). The 150 samples are displayed by small blue markers whilst the removed samples are displayed by round grey red crossed markers. In the bottom panel: $10^5$ MC validation runs. Blue markers show new validation samples successfully predicted by the bounds whilst red markers display samples which do not lay within the IPM bounds.</figcaption></figure>

<table><caption>Table 4: Area metric and Monte Carlo (MC) estimate of the true probability of error of IPMs optimized with program (8) and $N = 2000$ samples. The IPMs bounds are defined by combination of radial basis $n_u = n_l = 20$. We also present the a-posteriori and a-priori Scenario bounds from Theorems 3 and 2, computed for a confidence $\beta = 10^{-6}$.</caption><tbody><tr><td>N<sub>o</sub></td><td>0</td><td>10</td><td>50</td><td>100</td><td>125</td><td>150</td><td>175</td><td>200</td><td>225</td></tr><tr><td>CPU Time [s]</td><td>1.19</td><td>1.090</td><td>0.87</td><td>0.90</td><td>0.91</td><td>98.40</td><td>90.99</td><td>1.05</td><td>85.08</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Accuracy</td></tr><tr><td>A(θ<sup>*</sup>)</td><td>0.56</td><td>0.52</td><td>0.45</td><td>0.43</td><td>0.40</td><td>0.38</td><td>0.39</td><td>0.39</td><td>0.40</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Epistemic bounds on P<sub>f</sub></td></tr><tr><td>ε̄(n<sub>θ</sub>, N, β, N<sub>o</sub>) a-prio.</td><td>0.04</td><td>0.07</td><td>0.12</td><td>0.172</td><td>0.195</td><td>0.21</td><td>0.24</td><td>-</td><td>-</td></tr><tr><td>s<sub>N</sub><sup>*</sup></td><td>20</td><td>32</td><td>70</td><td>122</td><td>146</td><td>176</td><td>196</td><td>224</td><td>251</td></tr><tr><td>ε̄(s<sub>N</sub><sup>*</sup> + N<sub>o</sub>, N, β) a-post.</td><td>0.028</td><td>0.037</td><td>0.065</td><td>0.096</td><td>0.110</td><td>0.128</td><td>0.139</td><td>0.156</td><td>0.17</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Validation</td></tr><tr><td>P̂<sub>f</sub>(θ<sup>*</sup>) MC</td><td>0.0127</td><td>0.018</td><td>0.041</td><td>0.060</td><td>0.077</td><td>0.092</td><td>0.097</td><td>0.114</td><td>0.123</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 4: In the panel on the left: the optimized regression function $f(x; \theta^*)$ and two minimax layer IPMs constructed for the empirical levels $\gamma_1^*$ and $\gamma_{10}^*$. In the panel on the right: the probability box bounding the reliability of minimax layer IPMs for any $\gamma_i^*, i = 1, \dots, N$. The P-box holds with a confidence $\beta = 1 - 2 \cdot 10^{-4}$.</figcaption></figure>

### 5.1.3. Minimax layer

The Minimax layer program (10) is solved using a data set of $N = 2000$ samples and a regression function $f(x; \theta)$ defined by radial basis functions with scale parameter $\nu = 1$ and $n_\theta = 20$ fitting coefficients. Figure 4 presents the results. The optimized $f(x; \theta^*)$ is displayed by the red curve on the left panel and it has been used to construct two IPMsss $[f(x; \theta^*) - \gamma_1^*, f(x; \theta^*) + \gamma_1^*]$, and $[f(x; \theta^*) - \gamma_{10}^*, f(x; \theta^*) + \gamma_{10}^*]$, for the empirical levels $\gamma_1^*$ and $\gamma_{10}^*$, respectively. The empirical level $\gamma_1^* = \max_{i=1,...,N} |y_i - f(x; \theta^*)|$, corresponds to the value of the objective function at the optimum, i.e., is the maximum absolute vertical distance between the regression function and the samples. The level $\gamma_{10}^*$ corresponds to the $10^{th}$ highest difference $|y_i - f(x; \theta^*)|$. The reliability of the resulting IPMs is presented on the right panel by a reliability box constructed as described by Theorem 5.

The modeler can make use of the probability box to select the best IPM which compromises between accuracy and reliability. As an example consider the IPM for the level $\gamma_1^*$ (solid lines in the left panel) and $\gamma_{10}^*$ (dotted lines in the left panel). The accuracy of the predictor $[f(x; \theta^*) \pm \gamma_1^*]$ is given by an area $2\gamma_1^*$ which results in a larger area when compared to the area of the second predictor $2\gamma_{10}^*$. Theorem 5 guarantees that the reliability of the predictor $[f(x; \theta^*) \pm \gamma_1^*]$ is at least $R(\theta^*) > 0.99$. Although much more accurate, the reliability for a level $\gamma_{10}^*$ can not be better than 0.8152, but it is at least greater than 0.7206. A Monte Carlo estimation of the true reliability for each empirical level shows the validity of the derived bounds (see marked black CDF on the right panel).

Program (10) is easy to implement and if compared to other IPM optimization programs it has a small computational cost. This is due to the $N$ scenario constraints which are all compressed in a single objective function for the Minimax approach. Furthermore, powerful bounds on the cumulative distribution of the reliability of all the Minimax layer IPMs are obtained from a single run of the optimization program. One of the drawbacks of this implementation is that the optimized $f(x; \theta^*)$ may result in a poor fitting of Heteroscedastic data sets. In the proposed example, for instance, the optimization program is driven by the highest variance of $y$ in the left-hand side of the graph. This led to a good fit in that region and a relatively poor fit on the low-variance region.

### 5.1.4. Results for the proposed soft-constrained program

Figure 5 presents two examples of IPMs trained using program (34). Radial basis functions $\psi(x)$ and a degree $n_l = n_u = 10$ are selected for the bounds. In the top panel an IPM obtained for a cost of violations $\rho = 0.5$ is presented whilst the bottom panel presents an IPM trained with softer-constraints $\rho = 0.05$. Blue markers shows the $N = 1000$ training scenarios whilst the red x markers presents the scenarios which failed the validation and fall outside from the predictive intervals. It can be noticed that a higher cost of violations $\rho$ worsen the accuracy of the IPM leading predictive intervals which are (in average) wider. Specifically, a $\rho = 0.5$ leads to a $\mathcal{A} = 0.582$ whilst
<figure><img src="image_5.png" /><figcaption>Figure 5: Radial basis IPM bounds with $n_l = n_u = 10$ fitting coefficients optimized with program (34) for $N = 10^3$ samples (blue round markers) of DGM-1 having $\sigma = 3$. The red x markers show the failed validation samples. The top and bottom panels display the optimized IPM resulting from cost parameters $\rho = 0.5$ and $\rho = 0.05$, respectively.</figcaption></figure>

the area for a $\rho = 0.05$ is only $\mathcal{A} = 0.439$.

The IPM trained with hard constraints is more reliable. Specifically, we can guarantee that the reliability of the IPM trained for a $\rho = 0.5$ and $\rho = 0.05$ lay, respectively, in $R \in [0.96, 1.0]$ and $R \in [0.916, 0.98]$. In words, the model will predict correctly between 96 % and 100 % of the future scenarios. Differently, the second model will predict correctly at least 91.6 % of the future samples but not more that 98 % of them. Note that the gap between the lower and the upper bound quantify the epistemic uncertainty affecting our predictor and that arises due to a lack of scenarios $N$ and complexity of the IPM design given by $s_N^*$. Monte Carlo estimator of the true IPMs' reliability result $\hat{R} = 0.992$ and $\hat{R} = 0.956$ for $\rho = 0.5$ and $\rho = 0.05$, respectively. For both cases, $\hat{R}(\theta^*) \in [1 - \bar{\epsilon}, 1 - \epsilon]$ as prescribed by Scenario theory.

The coefficients of soft-constrained radial basis bounds of degree $n_\theta = 20$ are optimized using program (34) and for 50 values of $\rho$ between 0.01 and 0.5. We examine 20 training sets $\mathcal{D}_N$ of size $N \in [50, 10^3]$. Figures 6 and 7 present the results of this analysis. In Figure 6 the area $\mathcal{A}$, the scenario bounds [$\epsilon, \bar{\epsilon}$] from Theorem 4 and the MC estimate of $P_f$ are compared. Circle filled markers, blue circle markers and cross markers presents the result for increasing hardening of the scenario constraints (higher $\rho$ value). We can observe a trade off between area and both [$\epsilon, \bar{\epsilon}$] and $P_f$. Also, the failure probaiblity is positively correlated to $\epsilon$ and $\bar{\epsilon}$ although the correlation appears less marked with the upper bound. In Figure 7 we investigate further correlation between the estimated $P_f$ and the bounds. On the y-axis we present the normalized distance between the estimate value of $P_f$ and the lower bound given by $\frac{P_f - \epsilon}{\bar{\epsilon} - \epsilon}$. whilst on the x-axis we display the number of samples. Increasing the data set size $N$, the value of the true probability of failure gets closer to $\epsilon$ (blue markers). In contrast, smaller $N$ leads to a generalization interval which is tighter from above, i.e., $P_f$ results closer to $\bar{\epsilon}$ (larger red markers).

## 5.2. Case study II: identification of a black-box system-controller dynamic response model

For this case study, we seek an IPM characterizing the uncertain response of a system controlled by a non-collocated sensor/actuator presented in the NASA langley challenge on optimization under uncertainty [56]. The
<figure><img src="image_6.png" /><figcaption>Figure 6: The results for 1000 IPMs trained with program (34). The 3 marker types indicate different levels of relaxation for the constraints.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: The normalized difference between $(\frac{P_f - \epsilon}{\epsilon} - \epsilon)$ on y-axis against the 20 values of N (x-axis) and the 50 values of $\rho$ (iso-N markers). Interestingly, a higher N leads to a $P_f$ that is closer to the lower bound.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: An IPM for the dynamic response, $z_2(t)$, of the black-box model of a composite system controller with a non-collocated sensor/actuator pair [56].</figcaption></figure>

controller tries to suppress undesirable oscillations in the dynamic system which is affected by unknown sources of uncertainty and presented as a black-box model. We will focus our analysis on $z_2(t)$, which is responsible for the energy consumption of the controller [56]. The goal of case study II is to identify an IPM to model the uncertain response $z_2(t)$ from samples of it.

A data set of 100 responses is provided and each sample $z_2(t, i), t = 1, ..., T$ has a length $T = 5000$. This leads to a data set $\{\{z_2(t, i)\}_{i=1}^T\}_{t=1}^{100}$ and the total number of sample pairs is $5 \times 10^5$. We consider available a random subset of $N = 10^4$ samples and we use the remaining data points for validation of the model. Thee samples are used in program (34) to define $10^4$ constraints and identify an optimal IPM. A cost parameter $\rho = 10^5$ was selected and the bounds $f_u$ and $f_l$ defined by linear combinations of $n_u = n_l = 50$ radial basis functions having a scale parameter $\nu = 15$. Figure 8 presents the results. The IPM bounds are displayed by black solid lines, the red makers present the $10^4$ samples and the blue dotted lines show the 100 samples of $z_2$. The number of support scenarios is $s_{10^4}^* = 99$ and Theorem (4) ensures that the true IPM reliability is guaranteed to lay within the interval $[0.9852,0.9938]$ with a confidence $1 - 10^{-3}$. Sample-based estimation of the true reliability results in $\hat{R} = 0.9906$ and support the validity of the epistemic bounds.

Note that $\hat{R}$ is obtained averaging over the $N = 5 \times 10^5$ samples and not over the 100 sample responses. The IPM reliability indicates how likely a new sample point $z_2(t)$ will fall inside the bound prescribed by the IPM at time $t$. This probability differs from the probability

$$ \mathbb{P}[z_2(t) \in \mathcal{I}_N(t; \theta^*), \forall, t = 1, ..., T]. $$

that is the probability of $\mathcal{I}_N(t; \theta^*)$ containing a new realization of the process. We estimate this probability from the 100 samples of the process and only 74 out of the 100 sequences $z_2(1), ..., z_2(T)$ result fully contained within the bounds. A formal epistemic bound on this probability would be an important tool for model identification and reliable modelling of dynamic processes. Unfortunately, generalization bounds on the whole process can not be defined yet using Scenario theory. This is due to the non-stationary and dependent nature of the sequence and further research efforts are needed in this direction.

### 5.3. Case study III: frequency response of a car suspension arm damaged by a crack

For Case study III, we have to characterize the uncertain response of a numerical suspension arm crossed by a crack of unknown lengths and fixed position. First, we start defining a DGM for the process. A data set of 500 frequency response functions are collected from an extended FE model of the suspension arm by randomizing crack lengths in position one (see [57] for more details on the model). For reproducibility sake, we define the DGM of this
<figure><img src="image_9.png" /><figcaption>Figure 9: The first and second moments obtained from 500 frequency response functions of a cracked suspension arm.</figcaption></figure>

process by estimating the means and standard deviations of the frequency response over a set knots on the support $\mathcal{X}$. Figure 9 presents the mean values (black marked line) and mean ± the standard deviations (blue dashed and dotted lines). A data set $\mathcal{D}_N$ is generated by sampling $x_i \in \mathcal{X}$ uniformly over its support and then randomizing a frequency response $y_i$ based on mean and standard deviation and linear interpolation between the two closest knots. We assume samples $y$ to be Gaussian distributed. Note however that the knowledge of the underlying distribution is not used to construct an IPM and is reported in this section uniquely to help the reader reproduce the results.

### 5.3.1. Comparison between IPM programs and discussion

In this section, we compare the IPMs resulting from the data of case study III. As a figure of merit, we examine the accuracy of the predictors (given by their area metric) and the a-posteriori reliability bounds. We assume that $N = 2000$ scenarios have been collected from the unknown process and are available to characterize frequency response. Preliminary analysis shows that in this case study polynomial basis under-performs compared to radial basis functions. Also, a relatively high $n_\theta$ achieve better results. Thus, we define the bounds by linear combinations of $n_l = n_u = 80$ radial basis and fitting coefficients $\theta$. A scale parameter $\nu = 45$ is selected.

Table 5 presents the results of this analysis. The hard-constrained program (3), in the first row, led to the IPM with normalized area metric $\mathcal{A}(\theta^*) = 0.72$. The number of support constraints is $s_N^* = 160$. These are samples which define active constraints on the bounds and prevent the area from further improving. Theorem 2 is used to compute a-posteriori reliability bound guaranteeing a reliability of at least $R(\theta^*) > 0.881$ for this model and with a high confidence $1 - 10^{-6}$. In other words, the probability that the IPM will fail to correctly predict future frequency response values is not greater than 0.119. If the uncertainty affecting the reliability of the predictor is too large, the analyst is left with two options: 1) additional samples should be collected from the DGM (however this can be costly or not possible) and, 2) a simpler model with having a lower $n_\theta$, can be selected (however this likely leads to a lower accuracy).

If the average width of the IPM is too large, $\mathcal{A}$ can be reduced by discarding some of the samples. The second and third rows in Table 5 present the result of the IPM program (8) with $N_o = 200$ and $N_o = 400$ discarded constraints, respectively. As expected, the accuracy of the model can be substantially improved by removing samples and the normalized $\mathcal{A}$ result in 0.59 and 0.44 in the first and second case, respectively. However, the reliability guarantees we can prescribe for the IPM with discarded samples are much worse, resulting in only $R(\theta^*) > 1 - \epsilon(N, \beta, s_N^* + 200) = 0.724$ when 200 are discarded, and $R(\theta^*) > 1 - \epsilon(N, \beta, s_N^* + 400) = 0.519$ when 400 samples are removed. An alternative way to improve the IPM accuracy is to use our new soft-constrained IPM program (34). For a costs of violation $\rho = 0.05$ and $\rho = 0.01$, (34) results in an area of 0.62 and 0.52, respectively. Note that for the soft-constrained case, an upper bound on the reliability it is also available and when $\rho$ increases the bounds on the reliability get worse. This is due to the lower cost of constraint violations leading to larger $s_N^*$.
<figure><img src="image_10.png" /><figcaption>Figure 10: A comparison of radial basis IPMs trained using programs (3), discarding $N_o = 400$ samples, and the soft-constrained program with $\rho = 0.01$.</figcaption></figure>

<table><caption>Table 5: A comparison of computational time, accuracy and reliability bounds for the 4 IPMs trained using different programs and a data set $D_{2000}$ from the DGM of example III.</caption><thead><tr><th>Program Eq.</th><th>[min]</th><th>A(θ<sup>*</sup>)</th><th>R&#770;(θ<sup>*</sup>)</th><th>[1 - ε&#773;, 1 - ε&#773;]</th><th>Theor.</th></tr></thead><tbody><tr><td>(3)</td><td>3.4</td><td>0.72</td><td>0.924</td><td>[0.881, 1.0]</td><td>2</td></tr><tr><td>(8) N<sub>o</sub> = 200</td><td>23.4</td><td>0.59</td><td>0.859</td><td>[0.724, 1.0]</td><td>3</td></tr><tr><td>(8) N<sub>o</sub> = 400</td><td>19.3</td><td>0.44</td><td>0.731</td><td>[0.519, 1.0]</td><td>3</td></tr><tr><td>(10) γ<sub>20</sub><sup>*</sup></td><td>0.1</td><td>3.30</td><td>0.989</td><td>[0.973, 1.0]</td><td>5</td></tr><tr><td>(34) ρ = 0.05</td><td>12.0</td><td>0.62</td><td>0.874</td><td>[0.824, 0.91]</td><td>4</td></tr><tr><td>(34) ρ = 0.01</td><td>22.5</td><td>0.52</td><td>0.786</td><td>[0.726, 0.83]</td><td>4</td></tr></tbody></table>

The results for the Minimax layer IPM defined by the 20<sup>st</sup> highest empirical level $\gamma_{20}^*$ are presented on in the 4th row of Table 5. The area metric results $A = 3.3$, which is much larger compared to the other methods and that makes the resulting predictor very inaccurate. Nonetheless, the larger area leads to the model with the strongest reliability guarantees. The minimax methods is the most efficient and took only (approx) 6 seconds to complete the analysis. In contrast, the soft-constrained program and the hard-constrained program with discarded samples have the highest computational costs, greater than 10 minutes<fnref n="3" />. Figure 10 presents the bounds $I_N(x; \theta^*)$ optimized via program (3), program (8) optimally discarding $N_o = 400$ samples, and the soft-constrained program (34) with $\rho = 0.01$. The bounds are displayed by black solid lines, red dashed thinner lines and blue dashed line, respectively. The 2000 scenarios in the data set are presented by gray round markers.

## 6. Conclusion

Interval Predictor Models offer an interval-valued characterization of the uncertainty affecting processes while prescribing formal guarantees on the model's reliability for future, yet unobserved, data. In this paper, we introduced a new method for Interval Predictor Models identification which allows trading reliability, given by the probability that the interval will include future data, for model's accuracy, which is the average width of the predictive intervals. A comparison of epistemic reliability bounds, formally prescribed by Scenario optimization theory, is also proposed in this work. Scenario theory quantifies the lack of data uncertainty affecting the solutions of scenario programs and it was used in this work to equip the new IPM with epistemic upper and lower bounds on its reliability. Scenario-based bounds hold distribution-free, non-asymptotically and their width depends on the sample size available to characterize the uncertain process. Moreover, we show that simpler models generally lead to better bounds. This is due to the lower risk of over-fitting and uncertainty under-estimation associated with low-complexity designs. The proposed method has been numerically compared to traditional approaches using three realistic case studies. Different degrees of a lack of samples and models complexity are investigated, showing that the new method finds

<footnote n="3">³Note that the computational cost can be reduced of an order of magnitude by adopting linear programming solvers</footnote>
solutions with a good compromise between accuracy and reliability. Potential improvements for future research on
IPMs include:

Study the relationship between IPMs and different non-deterministic predictors such as e.g., Gaussian Process regressors, Fuzzy regressors, etc. The interested reader is reminded to the work of M. Feas, J. Sadeghi et al. [95] for preliminary comparison between Kiring models and IPMs.

Extend the scope of scenario theory to prescribe guarantees on the probability of an IPM containing a whole
realization for a non-stationary and time-correlated stochastic processes.

Improve the efficiency of the soft-constrained methods when a large number of samples are available. For
instance by large-scale linear programming strategies [96, 97].

Investigate non-convex IPM optimization programs and new epistemic bounds for non-convex cases. In par-
ticular, investigate bounds for non-convex problems with discarded samples and soft-constrained non-convex
programs.

Interval predictors models prove to be a powerful and versatile tool to characterize the uncertainty of a random
process and Scenario theory allowed the epistemic uncertainty affecting the reliability of the predictions to be for-
mally quantified. As a final concluding remark, note that scenario bounds apply in full generality to a wide range of
optimization problems and are therefore not limited to Interval Predictor Models. We hope that the general overview
of the theory, the scenario programs, and the epistemic bounds presented in this work will benefit the readers and
provide them with tools for the formal quantification of lack of data uncertainty in their field of study.

Acknowledgment

The authors would like to acknowledge the gracious support of this work through the local authorities under
grant agreement “ITEA-2018-17030-DayTiMe” on Digital Life-cycle Twins for predictive maintenance.

Replication package:

The data and scripts for the numerical analysis are publicly available at the GitHub repository:

https://github.com/Roberock/ScenarioIPM (IPM Matlab Class)

References

[1] M. Campi, G. Calafiore, S. Garatti, New results on the identification of interval predictor models, IFAC Proceedings Volumes 38 (2005) 106–111. 16th IFAC World Congress.

[2] L. G. Crespo, D. P. Giesy, S. P. Kenny, Interval predictor models with a formal characterization of uncertainty and reliability, in: 53rd IEEE Conference on Decision and Control, pp. 5991–5996.

[3] M. Eldred, L. Swiler, G. Tang, Mixed aleatory-epistemic uncertainty quantification with stochastic expansions and optimization-based interval estimation, Reliability Engineering & System Safety 96 (2011) 1092–1113. Quantification of Margins and Uncertainties.

[4] C. Wang, H. G. Matthies, A comparative study of two interval-random models for hybrid uncertainty propagation analysis, Mechanical Systems and Signal Processing 136 (2020) 106531.

[5] D. Wu, W. Gao, Hybrid uncertain static analysis with random and interval fields, Computer Methods in Applied Mechanics and Engineering 315 (2017) 222–246.

[6] R. Rocchetta, M. Broggi, E. Patelli, Do we have enough data? robust reliability via uncertainty quantification, Applied Mathematical Modelling 54 (2018) 710–721.

[7] R. Rocchetta, E. Zio, E. Patelli, A power-flow emulator approach for resilience assessment of repairable power grids subject to weather-induced failures and data deficiency, Applied Energy 210 (2018) 339–350.

[8] G. Shafer, A mathematical theory of evidence turns 40, International Journal of Approximate Reasoning 79 (2016) 7 – 25. 40 years of Research on Dempster-Shafer Theory.

[9] S. Ferson, V. Kreinovich, L. Ginzburg, D. Myers, K. Sentz, Constructing probability boxes and dempster-shafer structures, Technical Report Sandia National Labs (2003).

[10] D. Dubois, H. Prade, Possibility theory, probability theory and multiple-valued logics: A clarification, Annals of Mathematics and Artificial Intelligence 32 (2001) 35–66.

[11] P. Walley, Statistical Reasoning with Imprecise Probabilities, Chapman & Hall/CRC Monographs on Statistics & Applied Probability, Taylor & Francis, 1991.

[12] L. Zadeh, Fuzzy sets, Information and Control 8 (1965) 338 – 353.

[13] R. Boukezzoula, D. Coquin, A decision-making computational methodology for a class of type-2 fuzzy intervals: An interval-based approach, Information Sciences 510 (2020) 256 – 282.

[14] Z. Liu, Y. Liu, J. Dezert, F. Cuzzolin, Evidence combination based on credal belief redistribution for pattern classification, IEEE Transactions on Fuzzy Systems 28 (2020) 618–631.

[15] M. Beer, S. Ferson, V. Kreinovich, Imprecise probabilities in engineering analyses, Mechanical Systems and Signal Processing 37 (2013) 4–29.

[16] G. Calafiore, M. Campi, A learning theory approach to the construction of predictor models, Discrete and Continuous Dynamical Systems-Series A 9 (2002) 156–166.

[17] M. Campi, G. Calafiore, S. Garatti, Interval predictor models: Identification and reliability, Automatica 45 (2009) 382 – 392.
[18] V. N. Vapnik, Statistical Learning Theory, Wiley-Interscience, 1998.

[19] T. Alamo, R. Tempo, E. F., Statistical learning theory: A pack-based strategy for uncertain feasibility and optimization problems, in: V. D. Blondel, S. P. Boyd, H. Kimura (Eds.), Recent Advances in Learning and Control, Springer London, London, 2008, pp. 1-14.

[20] K. Margellos, M. Prandini, J. Lygeros, A compression learning perspective to scenario based optimization, in: 53rd IEEE Conference on Decision and Control, pp. 5997-6002.

[21] L. Romao, K. Margellos, A. Papachristodoulou, On the exact feasibility of convex scenario programs with discarded constraints, 2020.

[22] G. Calafiore, M. Campi, Uncertain convex programs: randomized solutions and confidence levels, Mathematical Programming 102 (2005) 25-46.

[23] G. C. Calafiore, M. C. Campi, The scenario approach to robust control design, IEEE Transactions on Automatic Control 51 (2006) 742-753.

[24] M. Campi, S. Garatti, The exact feasibility of randomized solutions of uncertain convex programs, SIAM Journal on Optimization 19 (2008) 1211-1230.

[25] M. C. Campi, S. Garatti, Wait-and-judge scenario optimization, Mathematical Programming 167 (2018) 155-189.

[26] F. A. Ramponi, M. C. Campi, Expected shortfall: Heuristics and certificates, European Journal of Operational Research 267 (2018) 1003-1013.

[27] M. C. Campi, S. Garatti, F. A. Ramponi, Non-convex scenario optimization with application to system identification, in: 2015 54th IEEE Conference on Decision and Control (CDC), pp. 4023-4028.

[28] M. C. Campi, S. Garatti, F. A. Ramponi, A general scenario theory for nonconvex optimization and decision making, IEEE Transactions on Automatic Control 63 (2018) 4067-4078.

[29] S. Grammatico, X. Zhang, K. Margellos, P. Goulart, J. Lygeros, A scenario approach for non-convex control design, IEEE Transactions on Automatic Control 61 (2016) 334-345.

[30] Y. Yang, C. Sutanto, Chance-constrained optimization for nonconvex programs using scenario-based methods, ISA Transactions 90 (2019) 157-168.

[31] A. Carè, S. Garatti, M. C. Campi, Scenario min-max optimization and the risk of empirical costs, SIAM Journal on Optimization 25 (2015) 2061-2080.

[32] L. G. Crespo, B. K. Colbert, S. P. Kenny, D. P. Giesy, On the quantification of aleatory and epistemic uncertainty using sliced-normal distributions, Systems & Control Letters 134 (2019) 104560.

[33] M. M. Siraj, M. B. Saltik, P. M. V. den Hof, S. Grammatico, Scenario-based robust optimization of water flooding in oil reservoirs enjoys probabilistic guarantees, IFAC-PapersOnLine 51 (2018) 102-107. 3rd IFAC Workshop on Automatic Control in Offshore Oil and Gas Production OOGP 2018.

[34] A. Carè, S. Garatti, M. Campi, The wait-and-judge scenario approach applied to antenna array design, Computational Management Science (2019).

[35] R. Rocchetta, L. Crespo, S. Kenny, Solution of the benchmark control problem by scenario optimization, in: Dynamic Systems and Control Conference, volume 2. https://doi.org/10.1115/DSCC2019-8949.

[36] R. Rocchetta, L. G. Crespo, S. P. Kenny, A scenario optimization approach to reliability-based design, Reliability Engineering & System Safety 196 (2020) 106755.

[37] R. Rocchetta, L. G. Crespo, An empirical approach to reliability-based design using scenario optimization, in: Proceedings of the 30th European Safety and Reliability Conference and the 15th Probabilistic Safety Assessment and Management Conference.

[38] M. C. Campi, S. Garatti, Scenario optimization with relaxation: a new tool for design and application to machine learning problems, aixr, 2020.

[39] D. Paccagnan, M. C. Campi, The scenario approach meets uncertain game theory and variational inequalities, 2019 IEEE 58th Conference on Decision and Control (CDC) (2019).

[40] L. Bugliari Armenio, L. Fagiano, E. Terzi, M. Farina, R. Scattolini, Scenario optimization for optimal training of Echo State Networks, arXiv e-prints (2019) arXiv:1912.01693.

[41] S. Garatti, M. C. Campi, Learning for control: a bayesian scenario approach, in: 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 1772-1777.

[42] J. Sadeghi, M. de Angelis, E. Patelli, Efficient training of interval neural networks for imprecise training data, Neural Networks 118 (2019) 338-351.

[43] M. Lacerda, S. Martins, E. Nepomuceno, Structure selection based on interval predictor model for recovering static nonlinearities from chaotic data, IET Control Theory & Applications 12 (2018).

[44] X. Xiao, P. Wang, J. Wang, Newton method for interval predictor model with sphere parameter set, in: ICCCS.

[45] J. Sadeghi, M. de Angelis, E. Patelli, Frequentist history matching with interval predictor models, Applied Mathematical Modelling 61 (2018) 29-48.

[46] L. G. Crespo, S. P. Kenny, D. P. Giesy, R. B. Norman, S. Blattnig, Application of Interval Predictor Models to Space Radiation Shielding, pp. 1-23.

[47] F. Cannizzaro, N. Impollonia, S. Caddemi, I. Caliò, Explicit dynamic response of damaged beams with application to uncertain and identification problems, Journal of Sound and Vibration 487 (2020) 115608.

[48] M. L. Fravolini, G. Costante, T. Yucelen, M. R. Napolitano, Interval prediction models for data-driven design of aerial vehicle's robust adaptive controllers, Journal of Guidance, Control, and Dynamics 43 (2020) 1512-1528.

[49] J. hong Wang, Application of interval predictor model into robust model predictive control, International Journal of Systems Science 50 (2019) 1434-1449.

[50] J. Sadeghi, M. de Angelis, E. Patelli, Robust propagation of probability boxes by interval predictor models, Structural Safety 82 (2020) 101889.

[51] L. G. Crespo, S. P. Kenny, D. P. Giesy, Staircase predictor models for reliability and risk analysis, Structural Safety 75 (2018) 35-44.

[52] E. Patelli, S. Tolo, H. Williams, J. Sadeghi, R. Rocchetta, M. Angelis, M. Broggi, Opencossan 2.0: an efficient computational toolbox for risk, reliability and resilience analysis, in: Proceedings of the joint ICVRAM ISUMA UNCERTAINTIES conference, 2018.

[53] A. Carè, S. Garatti, M. Campi, Fast: an algorithm for the scenario approach with reduced sample complexity, IFAC Proceedings Volumes 44 (2011) 9236-9241. 18th IFAC World Congress.

[54] S. Garatti, M. C. Campi, Risk and complexity in scenario optimization, Mathematical Programming https://doi.org/10.1007/s10107-019-01446-4 (2019).

[55] M. J. Lacerda, L. G. Crespo, Interval predictor models for data with measurement uncertainty, in: 2017 American Control Conference (ACC), pp. 1487-1492.

[56] L. G. Crespo, S. P. Kenny, The nasa langley challenge on optimization under uncertainty, Mechanical Systems and Signal Processing 152 (2021) 107405.
[57] R. Rocchetta, M. Broggi, Q. Huchet, E. Patelli, On-line bayesian model updating for structural health monitoring, Mechanical Systems and Signal Processing 103 (2018) 174–195.

[58] D. Kong, Y. Chen, N. Li, Gaussian process regression for tool wear prediction, Mechanical Systems and Signal Processing 104 (2018) 556–574.

[59] S. D. Veiga, A. Marrel, Gaussian process regression with linear inequality constraints, Reliability Engineering & System Safety 195 (2020) 106732.

[60] G. Falsone, Stochastic differential calculus for gaussian and non-gaussian noises: A critical review, Communications in Nonlinear Science and Numerical Simulation 56 (2018) 198–216.

[61] A. Jain, T. Nghiem, M. Morari, R. Mangharam, Learning and control using gaussian processes, in: 2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS), pp. 140–149.

[62] A. Garbuno-Inigo, F. A. DiazDelaO, K. M. Zuev, Gaussian process hyper-parameter estimation using parallel asymptotically independent markov sampling, 2015.

[63] Z. Zhang, X. Si, C. Hu, Y. Lei, Degradation data analysis and remaining useful life estimation: A review on wiener-process-based methods, European Journal of Operational Research 271 (2018) 775–796.

[64] R. Boukezzoula, D. Coquin, K. Jacq, A possibilistic regression based on gradual interval b-splines: Application for hyperspectral imaging lake sediments, Information Sciences 510 (2020) 183–204.

[65] C.-C. Chuang, Extended support vector interval regression networks for interval input-output data, Information Sciences 178 (2008) 871–891. Including Special Issue "Ambient Intelligence".

[66] L. V. Utkin, An imprecise extension of svm-based machine learning models, Neurocomputing 331 (2019) 18–32.

[67] Z. Yang, D. K. Lin, A. Zhang, Interval-valued data prediction via regularized artificial neural network, Neurocomputing 331 (2019) 336–345.

[68] M. Pratama, E. Lughofer, M. J. Er, S. Anavatti, C.-P. Lim, Data driven modelling based on recurrent interval-valued metacognitive scaffolding fuzzy neural network, Neurocomputing 262 (2017) 4–27. Online Real-Time Learning Strategies for Data Streams.

[69] M. M. Dannert, M. G. Faes, R. M. Fleury, A. Fau, U. Nackenhorst, D. Moens, Imprecise random field analysis for non-linear concrete damage analysis, Mechanical Systems and Signal Processing 150 (2021) 107343.

[70] M. Alamaniotis, A. Grelle, L. H. Tsoukalas, Regression to fuzziness method for estimation of remaining useful life in power plant components, Mechanical Systems and Signal Processing 48 (2014) 188–198.

[71] S.-E. Fang, Q.-H. Zhang, W.-X. Ren, An interval model updating strategy using interval response surface models, Mechanical Systems and Signal Processing 60-61 (2015) 909–927.

[72] L. G. Crespo, S. P. Kenny, D. P. Giesy, Random predictor models for rigorous uncertainty quantification, International Journal for Uncertainty Quantification 5 (2015) 469–489.

[73] L. Wang, J. M. Mendel, Fuzzy basis functions, universal approximation, and orthogonal least-squares learning, IEEE Transactions on Neural Networks 3 (1992) 807–814.

[74] A. Graps, An introduction to wavelets, IEEE Comput. Sci. Eng. 2 (1995) 50–61.

[75] S. Garatti, M. Campi, A. Carè, On a class of interval predictor models with universal reliability, Automatica 110 (2019) 108542.

[76] L. G. Crespo, S. P. Kenny, D. P. Giesy, Interval Predictor Models With a Linear Parameter Dependency, Journal of Verification, Validation and Uncertainty Quantification 1 (2016). 021007.

[77] L. Cizelj, B. Mavko, H. Riesch-Oppermann, Application of first and second order reliability methods in the safety assessment of cracked steam generator tubing, Nuclear Engineering and Design (1994) 359–368.

[78] H.-S. Li, S.-K. Au, Design optimization using subset simulation algorithm, Structural Safety 32 (2010) 384–392. Modeling and Analysis of Rare and Imprecise Information.

[79] D. Altieri, E. Tubaldi, M. De Angelis, E. Patelli, A. Dall'Asta, Reliability-based optimal design of nonlinear viscous dampers for the seismic protection of structural systems, Bulletin of Earthquake Engineering 16 (2018) 963–982.

[80] M. de Angelis, E. Patelli, M. Beer, Advanced line sampling for efficient robust reliability analysis, Structural Safety 52 (2015) 170–182. Engineering Analyses with Vague and Imprecise Information.

[81] J. C. Medina, A. A. Taflanidis, Adaptive importance sampling for optimization under uncertainty problems, Computer Methods in Applied Mechanics and Engineering 279 (2014) 133–162.

[82] F. Grooteman, An adaptive directional importance sampling method for structural reliability, Probabilistic Engineering Mechanics 26 (2011) 134–141.

[83] K. Deb, S. Gupta, D. Daum, J. Branke, A. K. Mall, D. Padmanabhan, Reliability-based optimization using evolutionary algorithms, IEEE Transactions on Evolutionary Computation 13 (2009) 1054–1074.

[84] Z. Meng, H. Zhou, G. Li, D. Yang, A decoupled approach for non-probabilistic reliability-based design optimization, Computers & Structures 175 (2016) 65–73.

[85] Y. Kanno, A data-driven approach to non-parametric reliability-based design optimization of structures with uncertain load, Structural and Multidisciplinary Optimization (2019).

[86] L. Jaulin, E. Walter, Guaranteed robust nonlinear minimax estimation, IEEE Transactions on Automatic Control 47 (2002) 1857–1864.

[87] M. C. Campi, S. Garatti, A sampling-and-discarding approach to chance-constrained optimization: Feasibility and optimality, Journal of Optimization Theory and Applications 148 (2011) 257–280.

[88] R. Schöbi, B. Sudret, Global sensitivity analysis in the context of imprecise probabilities (p-boxes) using sparse polynomial chaos expansions, Reliability Engineering & System Safety 187 (2019) 129–141. Sensitivity Analysis of Model Output.

[89] W. Zhu, N. Chen, J. Liu, M. Beer, A probability-box-based method for propagation of multiple types of epistemic uncertainties and its application on composite structural-acoustic system, Mechanical Systems and Signal Processing 149 (2021) 107184.

[90] R. Rocchetta, E. Patelli, Assessment of power grid vulnerabilities accounting for stochastic loads and model imprecision, International Journal of Electrical Power & Energy Systems 98 (2018) 219–232.

[91] R. Rocchetta, E. Patelli, A post-contingency power flow emulator for generalized probabilistic risks assessment of power grids, Reliability Engineering & System Safety 197 (2020) 106817.

[92] S. Garatti, M. C. Campi, Complexity-based modulation of the data-set in scenario optimization, in: 2019 18th European Control Conference (ECC), pp. 1386–1391.

[93] X. Geng, L. Xie, Data-driven decision making with probabilistic guarantees (part 1): A schematic overview of chance-constrained optimization, arXiv:1903.10621, 2019.

[94] M. A. P. K. D. Chatterjee, R. Banavar, Scenario approach for minmax optimization with emphasis on the nonconvex case: positive results and caveats, arXiv:1906.01476, 2019.

[95] M. Faes, J. Sadeghi, M. Broggi, M. De Angelis, E. Patelli, M. Beer, D. Moens, On the robust estimation of small failure probabilities for strong nonlinear models 1 (2019).
[96] Y. Nesterov, Subgradient methods for huge-scale optimization problems, Mathematical Programming 146 (2013).

[97] G. Lancia, P. Serafini, Large-Scale Linear Programming, pp. 67–74.