Roberto Rocchetta
email: r.rocchetta@tue.nl

Qi Gao
email: q.gao@philips.com

Milan Petkovic
email: milan.petkovic@philips.com

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

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

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.

Introduction

An Interval Predictor Model is a rule that assigns to explanatory inputs ranges of plausible outcomes of a dependent variable [START_REF] Campi | New results on the identification of interval predictor models[END_REF]. The uncertainty in the process is quantified via non-probabilistic, interval-valued, predictions [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF]. Uncertainties can be modelled probabilistically, non-probabilistically or by a mixture of these two approaches [START_REF] Eldred | Mixed aleatory-epistemic uncertainty quantification with stochastic expansions and optimization-based interval estimation[END_REF][START_REF] Wang | A comparative study of two interval-random models for hybrid uncertainty propagation analysis[END_REF][START_REF] Wu | Hybrid uncertain static analysis with random and interval fields[END_REF]. 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 [START_REF] Rocchetta | Do we have enough data? robust reliability via uncertainty quantification[END_REF]. 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 [START_REF] Rocchetta | A power-flow emulator approach for resilience assessment of repairable power grids subject to weatherinduced failures and data deficiency[END_REF]. Evidence theory [START_REF] Shafer | A mathematical theory of evidence turns 40[END_REF][START_REF] Ferson | Constructing probability boxes and dempster-shafer structures[END_REF], Possibility theory [START_REF] Dubois | Possibility theory, probability theory and multiple-valued logics: A clarification[END_REF], Credal sets, Fuzzy sets and Ambiguity sets theory [START_REF] Walley | Statistical Reasoning with Imprecise Probabilities[END_REF][START_REF] Zadeh | Fuzzy sets[END_REF][START_REF] Boukezzoula | A decision-making computational methodology for a class of type-2 fuzzy intervals: An interval-based approach[END_REF][START_REF] Liu | Evidence combination based on credal belief redistribution for pattern classification[END_REF], are some of the most commonly used paradigms for this [START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF].

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 [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF][START_REF] Campi | Interval predictor models: Identification and reliability[END_REF]. 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 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 [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Alamo | Statistical learning theory: A pack-based strategy for uncertain feasibility and optimization problems[END_REF] and compression learning [START_REF] Margellos | A compression learning perspective to scenario based optimization[END_REF][START_REF] Romao | On the exact feasibility of convex scenario programs with discarded constraints[END_REF], 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, [START_REF] Calafiore | Uncertain convex programs: randomized solutions and confidence levels[END_REF][START_REF] Calafiore | The scenario approach to robust control design[END_REF][START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF][START_REF] Campi | Wait-and-judge scenario optimization[END_REF][START_REF] Ramponi | Expected shortfall: Heuristics and certificates[END_REF], and was recently extended to non-convex and soft-constrained cases, [START_REF] Campi | Non-convex scenario optimization with application to system identification[END_REF][START_REF] Campi | A general scenario theory for nonconvex optimization and decision making[END_REF][START_REF] Grammatico | A scenario approach for non-convex control design[END_REF][START_REF] Yang | Chance-constrained optimization for nonconvex programs using scenario-based methods[END_REF]. For the solution of a convex problem, Campi and Garatti [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF] 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 [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF], 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. [START_REF] Carè | Scenario min-max optimization and the risk of empirical costs[END_REF] 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 [START_REF] Crespo | On the quantification of aleatory and epistemic uncertainty using sliced-normal distributions[END_REF], a water flooding optimization [START_REF] Siraj | Scenario-based robust optimization of water flooding in oil reservoirs enjoys probabilistic guarantees[END_REF], an antenna arrays design [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF], optimal control and reliability-based design-optimization problems [START_REF] Rocchetta | Solution of the benchmark control problem by scenario optimization[END_REF][START_REF] Rocchetta | A scenario optimization approach to reliability-based design[END_REF][START_REF] Rocchetta | An empirical approach to reliability-based design using scenario optimization[END_REF]. More recently, the generalization of Support Vector Machine models [START_REF] Campi | Scenario optimization with relaxation: a new tool for design and application to machine learning problems[END_REF], Game-theoretic policies [START_REF] Paccagnan | The scenario approach meets uncertain game theory and variational inequalities[END_REF] and other machine learning tools [START_REF] Bugliari Armenio | Scenario optimization for optimal training of Echo State Networks[END_REF][START_REF] Garatti | Learning for control: a bayesian scenario approach[END_REF][START_REF] Sadeghi | Efficient training of interval neural networks for imprecise training data[END_REF][START_REF] Lacerda | Structure selection based on interval predictor model for recovering static nonlinearities from chaotic data[END_REF][START_REF] Xiao | Newton method for interval predictor model with sphere parameter set[END_REF] 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 [START_REF] Sadeghi | Frequentist history matching with interval predictor models[END_REF], space radiation shielding analysis [START_REF] Crespo | Application of Interval Predictor Models to Space Radiation Shielding[END_REF], damaged beams dynamic response modeling [START_REF] Cannizzaro | Explicit dynamic response of damaged beams with application to uncertain and identification problems[END_REF], design of adaptive controller for Aerial Vehicle's [START_REF] Fravolini | Interval prediction models for data-driven design of aerial vehicle's robust adaptive controllers[END_REF], robust predictive control [START_REF] Wang | Application of interval predictor model into robust model predictive control[END_REF], robust uncertainty propagation [START_REF] Sadeghi | Robust propagation of probability boxes by interval predictor models[END_REF], and for the detection and removal outliers [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF]. In [START_REF] Crespo | Staircase predictor models for reliability and risk analysis[END_REF] a random predictor model was introduced and compared with IPMs. The interested reader is reminded to, e.g., [START_REF] Rocchetta | A scenario optimization approach to reliability-based design[END_REF][START_REF] Patelli | Opencossan 2.0: an efficient computational toolbox for risk, reliability and resilience analysis[END_REF][START_REF] Carè | Fast: an algorithm for the scenario approach with reduced sample complexity[END_REF] 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, [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Margellos | A compression learning perspective to scenario based optimization[END_REF][START_REF] Romao | On the exact feasibility of convex scenario programs with discarded constraints[END_REF][START_REF] Grammatico | A scenario approach for non-convex control design[END_REF]. 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 [START_REF] Garatti | Risk and complexity in scenario optimization[END_REF] 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 [START_REF] Lacerda | Interval predictor models for data with measurement uncertainty[END_REF]; ii) on the identification of a model for the dynamic response of a black-box system-controller affected by uncertainty [START_REF] Crespo | The nasa langley challenge on optimization under uncertainty[END_REF]; 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 [START_REF] Rocchetta | On-line bayesian model updating for structural health monitoring[END_REF].

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.

Interval Predictor Model

A stochastic data generating mechanism is defined by a probability space (∆, F, P), where ∆ is the event space equipped with a σ-algebra F and a probability measure P. The DGM generates samples (x(t), y(t)) from an unknown random process y(t) = f (x(t)) + E, where E is a random noise term and a pair δ(t) = (x(t), y(t)) ∈ ∆ is composed 2 by a vector of explanatory variables x(t) ∈ X ⊆ R n x and a vector of dependent variables y(t) ∈ Y ⊆ R n y such that ∆ = X × 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 D N = {x i , y i } N i=1 available for the characterization of the process. The input-output samples δ i = (x i , y i ) ∈ 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 f (x) : R n x → R n y , approximating the true f (x). Although 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 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 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 [START_REF] Kong | Gaussian process regression for tool wear prediction[END_REF][START_REF] Veiga | Gaussian process regression with linear inequality constraints[END_REF][START_REF] Falsone | Stochastic differential calculus for gaussian and non-gaussian noises: A critical review[END_REF][START_REF] Jain | Learning and control using gaussian processes[END_REF][START_REF] Garbuno-Inigo | Gaussian process hyper-parameter estimation using parallel asymptotically independent markov sampling[END_REF], Wiener and Poisson processes [START_REF] Zhang | Degradation data analysis and remaining useful life estimation: A review on wiener-process-based methods[END_REF], Interval-valued B-splines [START_REF] Boukezzoula | A possibilistic regression based on gradual interval b-splines: Application for hyperspectral imaging lake sediments[END_REF], support vector regression models [START_REF] Chuang | Extended support vector interval regression networks for interval input-output data[END_REF][START_REF] Utkin | An imprecise extension of svm-based machine learning models[END_REF], neural networks [START_REF] Yang | Interval-valued data prediction via regularized artificial neural network[END_REF][START_REF] Pratama | Data driven modelling based on recurrent interval-valued metacognitive scaffolding fuzzy neural network[END_REF], random fields [START_REF] Dannert | Imprecise random field analysis for non-linear concrete damage analysis[END_REF], Fuzzy regression models [START_REF] Alamaniotis | Regression to fuzziness method for estimation of remaining useful life in power plant components[END_REF], interval response surfaces [START_REF] Fang | An interval model updating strategy using interval response surface models[END_REF], or Random predictive models [START_REF] Crespo | Staircase predictor models for reliability and risk analysis[END_REF][START_REF] Crespo | Random predictor models for rigorous uncertainty quantification[END_REF].

Interval Predictive Model offers an interval-valued characterization of the uncertainty affecting a stochastic data generating process [START_REF] Crespo | Application of Interval Predictor Models to Space Radiation Shielding[END_REF] by defining a rule I N (x; θ), which assigns to a vector explanatory variable x a lower and upper bound I N (x; θ) = [ f l (x; θ), f u (x; θ)]. A vector of fitting coefficients θ defines two functions f l (x; θ) and f u (x; θ) bounding the dependent quantity of interest y. The resulting interval I N (x; θ), 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 ∈ I N (x; θ). 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 δ = 2. This has not to be considered as a limiting assumption. In fact, the dimension n δ 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.

Stationary and independent process

For simplicity sake, we assume the process δ 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 δ (δ 1+τ , ..., δ t+τ , ...) = F δ (δ 1 , ..., δ t , ...) ∀ τ, t ∈ R.

In other words, a stochastic process is stationary if the distribution of δ(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 δ t , t = 1, 2, ... rules out the presence of inter-time correlations:

F δ 1 ,...δ t (δ 1 , ..., δ t ) = F δ 1 (δ 1 ) • ... • F δ t (δ t ) ∀ δ i .
Independence in x(t) is an assumption which is generally met in many identification problems such as in econometric, pattern recognition, learning theory where no dynamics is present [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF]. 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, [START_REF] Campi | Interval predictor models: Identification and reliability[END_REF] 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.

Basis functions

The functions f u and f l are represented as a linear combination of fitting coefficients and basis function:

f u (x; θ) = n u j=0 u j ψ j (x), f l (x; θ) = n l j=0 l j ψ j (x), (1)
where θ = [u 0 , .., u n u , l 0 , ..., l n l ] ∈ R n θ is the vector of fitting coefficients, ψ 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 θ = n l + n u . Table 1 presents examples of basis function ψ j (x) which can be used to construct an IPM. Different basis, such as Fuzzy [START_REF] Wang | Fuzzy basis functions, universal approximation, and orthogonal least-squares learning[END_REF] or Wavelets basis [START_REF] Graps | An introduction to wavelets[END_REF], can be considered to improve the effectiveness of the model. The problem of selecting a good ψ j (x) is data and process dependent and in the remainder of this work ψ j (x) will be selected on an empirical basis. 3

Accuracy of an IPM

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

A(θ) = X f u (x; θ) -f l (x; θ)dx. (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 , ..., x m ∈ X in descending order and equally spaced by a distance ∆ x = x j+1x j . Then, the points are used for the estimation of the integral as follows,

A(θ) ≈ ∆ x m j=1 f u (x j-1 ; θ) -f u (x j ; θ) 2 -∆ x m j=1 f l (x j-1 ; θ) -f l (x j ; θ) 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 X is a one dimensional integral [x, x] and polynomial basis are selected, the area can be expressed analytically as follows:

A(θ) = n u j=0 u j x j+1 j x x - n l j=0 l j x j+1 j x x .
Similarly, if radial basis functions ψ j (x) = e -ν(x-c j ) 2 define the upper and lower bounds, the area metric is proportional to:

A(θ) ∝ n u j=0 u j er f √ ν(x -c j ) x x - n l j=0 l j er f √ ν(x -c j ) x x .
were er f (•) is the error function and c j are the knots values in X defining the centers of the radial basis ψ j (x). Note that the area A(θ) 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 f (x). This IPM is considered very precise because it yields prediction intervals that are very narrow, i.e., f (x) ± 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.

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 A(θ) while constraining all the available samples within I N (x; θ). Table 1: A list of commonly adopted basis functions. The parameter c i are user-defined knots values in X and ν is a scale parameter defining the spread of a radial basis from its center c i .

Name

Basis function

Polynomial basis ψ j (x i ) = x j i Radial (Gaussian) basis ψ j (x i ) = exp(-ν||c jx i ||) Fourier basis ψ j (x i ) = cos(c j x i ) for i odd Fourier basis ψ j (x i ) = sin(c j x i ) for i even Sigmoid/Soft Thresholds ψ j (x i ) = IPM program 1. (Full data enclosure with no exceptions [START_REF] Campi | Interval predictor models: Identification and reliability[END_REF])

θ = arg min θ A(θ) s.t. g(δ i ; θ) ≤ 0, i = 1, ..., N, (3)
where θ is the optimized vector of fitting coefficients (the IPM design), and the function g(δ; θ) defines joint constraints on the upper and lower bounds.

The random samples in the data set D N define N constraints in Eq. ( 3). The constraints are given by:

g(δ i ; θ) = max {g u (δ i ; θ), g l (δ i ; θ)} , i = 1, ..., N, (4)
where

g u (δ i ; θ) = y i -f u (x i ; θ), i = 1, ..., N, (5)
define the data-enclosure constraints for the upper bound f u , and

g l (δ i ; θ) = f l (x i ; θ) -y i , i = 1, ..., N, (6)
define the constraints for the lower bound. To ensure f u (x; θ) > f l (x; θ) for all x ∈ X, the following upper dominance constraint is included in the optimization program:

g d (δ; θ) = sup x∈X [ f l (x; θ) -f u (x; θ)] ≤ 0. ( 7
)
If ψ j (x) is a polynomial basis, a sum-of-squares formulation of this constraint can be used, see e.g. [START_REF] Lacerda | Interval predictor models for data with measurement uncertainty[END_REF] 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; θ)f u (x; θ), evaluated on an arbitrary number of points x 0 , x 1 , ..., x m ∈ X.

The IPM resulting from program (3) tightly encloses all the samples (x i , y i ) ∈ D N . This is due to the N constraints defined by ( 5) and ( 6), which must be satisfied with no exceptions for θ 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 D N . To this end, the authors of [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF] proposed to optimally discard samples from D N while minimizing A(θ) and enclosing between the bounds the remaining samples [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF]. The optimization program in Eq. ( 8) summarizes this approach.

IPM program 2. (Optimal removal of N o samples [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF])

θ = arg min θ A(θ) s.t. g(δ i ; θ) ≤ 0, ∀δ i ∈ D N-N o , (8)
where N o is the number of discarded samples, the set

D N-N o = δ i ∈ D N : Fg(δ i ;θ) ≤ N -N o N , (9)
contains the remaining N -N o samples defining constraints on f u and f l , and Fg(δ i ;θ) 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 [START_REF] Garatti | On a class of interval predictor models with universal reliability[END_REF] and named Minimax layer approach. This method, differently from (3) and ( 8), minimizes the maximum absolute distance between a regression function f (x i ; θ) and the samples y i . [START_REF] Garatti | On a class of interval predictor models with universal reliability[END_REF])

IPM program 3. (Minimax layer

θ = arg min θ max i∈{1,...,N} | f (x i ; θ) -y i |, (10)
where f (x i ; θ) = ψ T (x i )θ is a regression function defined by a linear combination of basis functions ψ T and fitting coefficients θ as presented in equations (1).

The Minimax layer IPM optimizes f (x i ; θ) by minimizing the L ∞ -norm be the errors [START_REF] Garatti | On a class of interval predictor models with universal reliability[END_REF], 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 ; θ ) is obtained, a constant width IPM is constructed as follows [START_REF] Garatti | On a class of interval predictor models with universal reliability[END_REF]:

I N (x; θ ) = [ f (x; θ ) -γ i ; f (x; θ ) + γ i ],
where

γ i = | f (x i ; θ ) -y i |
is the i th highest empirical cost value, i.e., the i th largest difference between the regression function f (x i ; θ ) and the samples of the dependent quantity. This IPM has a constant width of 2γ i and the highest empirical cost value γ 1 = max i∈{1,...,N}

| f (x i ; θ )y i | renders the IPM with the lowest accuracy (largest bounds) and highest reliability.

Convexity of the reviewed IPM programs

The functions f u and f l have a linear dependency with the optimization variables θ. The area metric also displays a linear dependency in θ. In fact, by assuming where f u (x; θ) ≥ f l (x; θ) for all x ∈ X:

A(θ) = X         n u j=0 u j ψ j (x) - n l j=0 l j ψ j (x)         dx = n u j=0 u j X ψ j (x)dx - n l j=0 l j X ψ j (x)dx.
The reviewed programs are convex with respect to the design parameters θ and can be solved very efficiently using a linear programming algorithm 1 . For further details on linear programming implementation of IPM see for instance [START_REF] Crespo | Interval Predictor Models With a Linear Parameter Dependency[END_REF].

Reliability of an IPM

The reliability of an IPM is the probability that I N (x; θ) will contain future samples from the DGM and is mathematically defined as follow:

R(θ) = P[(x, y) ∈ ∆ : y ∈ I N (x; θ)]. (11)
Similarly, the probability of error (or failure probability) of an IPM is given by:

P f (θ) = 1 -R(θ) = P[(x, y) ∈ ∆ : y I N (x; θ)]. (12)
Note that the true values of R(θ) and P f (θ) are unknown in practice. This is due to a limited availability of samples (x, y) ∈ D N and, consequently, an epistemic uncertainty affecting the probability P generating the data. Sample-based approximations of R(θ) and P f (θ) can be obtained from the available data as follows:

R(θ) = n i=1 1{y i ∈ I N (x i ; θ)} n , P f (θ) = n i=1 1{y i I N (x i ; θ)} n , ( 13
)
where n is the number of new samples used for the estimation and 1{c i } is the indicator function for a condition {c i }:

1{c i } =        1 if c i occurs 0 otherwise .
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→∞ P |R(θ) -R(θ)| ≥ = 0, ( 14
)
where 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 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

P[R(θ) -Rm (θ) ≥ ] = n(p-) k=0 n k p k (1 -p) n-k .
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(θ) and P f (θ) by prescribing a formal mathematical expressions for the bound . 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 θ of a scenario optimization problem. Scenario optimization programs and formal reliability bounds prescribed via Scenario theory will be presented next.

A summary of Scenario theoretic epistemic reliability bounds

Consider a chance-constrained reliability program:

θ = arg min θ∈Θ {J(θ) : R(θ) > α} , (15)
where θ ∈ Θ ⊆ R n θ is a vector of design variables constrained in a closed set Θ, J : Θ → R is a cost function and the reliability constraint R(θ) = P[δ ∈ ∆ : g(δ; θ) ≤ 0], enforces non-positiveness of g with a probability at least α ∈ [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 [START_REF] Cizelj | Application of first and second order reliability methods in the safety assessment of cracked steam generator tubing[END_REF], Subset Simulation [START_REF] Li | Design optimization using subset simulation algorithm[END_REF][START_REF] Altieri | Reliability-based optimal design of nonlinear viscous dampers for the seismic protection of structural systems[END_REF], Line Sampling [START_REF] De Angelis | Advanced line sampling for efficient robust reliability analysis[END_REF], Importance Sampling [START_REF] Medina | Adaptive importance sampling for optimization under uncertainty problems[END_REF][START_REF] Grooteman | An adaptive directional importance sampling method for structural reliability[END_REF] or other means [START_REF] Deb | Reliability-based optimization using evolutionary algorithms[END_REF][START_REF] Meng | A decoupled approach for non-probabilistic reliability-based design optimization[END_REF][START_REF] Kanno | A data-driven approach to non-parametric reliability-based design optimization of structures with uncertain load[END_REF]. The vast majority of the existing methods require a model of the uncertainty for δ to carry on with the analysis. However, this can lead to an overconfidence in the reliability R(θ ) and to an under-estimation of the true level of uncertainty affecting the performance of θ , especially under a lack of samples [START_REF] Rocchetta | Do we have enough data? robust reliability via uncertainty quantification[END_REF][START_REF] Rocchetta | Solution of the benchmark control problem by scenario optimization[END_REF].

Differently from traditional methods, Scenario programs gives a feasible solutions for Eq. ( 15) directly from N samples of δ and without assumptions on the family of P[•] generating the data. Scenario theory prescribes formal epistemic bounds on the reliability R(θ ) which quantify the worst-case reliability among the possible P[•] which are consistent with the observed samples. The bounds hold distribution-free, i.e., with no distributional assumptions on the underlying P[•], 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.

Scenario optimization programs

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

θ = arg min θ∈Θ {J(θ) : g(δ i ; θ) ≤ 0, i = 1, .., N, } (16)
where g(δ; θ) : Θ × ∆ → R is a constraint function defining N samples constraints. The constraints in [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF] are hard because an optimal design θ to be feasible must satisfy g(δ i ; θ ) ≤ 0, for all i = 1, ..., 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 θ and define, respectively, the bounds-enclosure constraints and the area metric objective function.

In may cases, ( 16) can lead to a feasibility issue if θ ∈ Θ : g(δ; θ) ≤ 0, for all δ ∈ D N . To amend for this deficiency, a scenario program with discarded samples can be considered as presented by [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF]. Scenario program 2. (Hard-constrained, discarded samples) A feasible sampling-and-discarding reformulation of program ( 16) is given by:

θ = arg min θ∈Θ {J(θ) : g(δ i ; θ) ≤ 0, δ i ∈ D N \ O} , (17)
where D N \ O is the set difference between the data and a set of discarded samples O. The cardinality of the set of discarded scenarios is N o = |O|.

The set 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 [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF], enjoys an improved cost J(θ ). 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 D N \ 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 [START_REF] Garatti | Risk and complexity in scenario optimization[END_REF].

Scenario program 3. (Soft-constrained) An feasible reformulation of ( 16) is given by:

θ , ζ = arg min θ∈Θ ζ∈R N +        J(θ) + ρ N i=1 ζ i : g(δ i ; θ) ≤ ζ i , δ i ∈ D N        , ( 18
)
where ζ is a N-dimensional vector of non-negative slack variables and ρ defines the cost of violations. A ζ i = 0 means that the hard-constraint imposed by the i th sample is satisfied, i.e., g(δ i ; θ) ≤ 0. On the other hand, a ζ 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 ρ → ∞ 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 epigraphyc formulation of the problem is defined as follows:

Scenario program 4. (minimax or L ∞ minimization) θ , γ = arg min θ∈Θ γ∈R γ : max i∈{1,...,N} g(δ i ; θ) ≤ γ , ( 19
)
where γ is the minimized L ∞ norm (maximum) of g(δ i ; θ ), i = 1, ..., N at the optimum θ .

The programs seek a design θ 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(δ i ; θ) = | f (x i ; θ)y i | and is, therefore, convex in the design variables θ. For earlier works on minimax robust approaches see, e.g., [START_REF] Jaulin | Guaranteed robust nonlinear minimax estimation[END_REF] and [START_REF] Wang | Application of interval predictor model into robust model predictive control[END_REF].

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

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

P f (θ ) = P[δ ∈ ∆ : θ Θ δ ], (20)
is called violation probability, or error probability. The quantity Θ δ is the feasibility set induced by δ. Given a reliability parameter ∈ [0, 1] a solution θ is called -robust if P f (θ ) < .

The quantity P f (θ ) gives the 'exact' probability that a solution θ will fail to meet the desired requirements, i.e., the probability that, given a new observation δ ∈ ∆, the constraint imposed by g(θ , δ) in correspondence of θ will not be satisfied. The feasibility set induced by the constraints in program (1) is defined as follows:

Θ δ = {θ ∈ Θ : g(δ; θ) ≤ 0} . (21)
Similarly, the feasibility set induced by the constraint in program ( 4) is given by:

Θ δ = θ ∈ Θ : g(δ; θ) ≤ γ . (22)
Note that for IPM programs the violation probability corresponds to the probability of the interval failing to predict a future observation correctly. Existence can be enforced by constraining the optimization on a compact domain Θ, this avoids problems which may arise when the objective function decreases when θ drifts to infinity. Concerning the uniqueness of the solution, if multiple optimal solutions exist in Θ a tie-break rule can be implemented, e.g. selecting the solution with minimum norm or optimizing additional convex functions in θ.

Assumption 2. (Non-degeneracy) For any positive integer N ∈ N 0 and scenario set 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(δ i , θ ) < 0, yield a new optimal solution having a smaller cost J(θ ) [START_REF] Rocchetta | A scenario optimization approach to reliability-based design[END_REF].

Assumption 3. (Non-accumulation) For every θ ∈ Θ the statement P[δ : g(δ; θ) = a] = 0 for any a ∈ R holds true.

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

An overview of epistemic reliability bounds

Scenario theory gives formal epistemic bounds on the probability of an optimized scenario design θ failing to comply with new observations, i.e., formal bounds on P f (θ ) and R(θ ). These bounds are also know 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.

Bounds for hard-constrained convex problems

Theorem 1. [24, Theorems 2.4] Consider a convex scenario program defined as in [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF]. Under assumptions 1, 2, stationary P and IID samples in D N and N > n θ the distribution of P f (θ ) is bounded by a Beta distribution:

P N [P f (θ ) ≤ ] ≥ 1 - n θ -1 k=0 N k k (1 -) N-k = 1 -β, (23)
where n θ is the number of design variables and β ∈ [0, 1] is a confidence parameter.

The parameter gives an upper bound on the probability of error and selecting a confidence parameter β, it can be expressed as a function (N, β, n θ ) of the number of design variables n θ , the number of samples N, and a confidence level β 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 θ of Eq. ( 16), i.e., it is a-priori computed. Theorem 1 was derived from the Helly's theorem and showing that s N ≤ n θ for hard-constrained convex scenario program. In other words, the complexity of the solution of a convex scenario program is always capped by n θ and the solution of the problem always admits a compressed representation of dimension n θ . For optimization program like in Eq. ( 16), when fullysupported, we have that s N = n θ 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 θ can be often reconstructed by a number of constraints s N < n θ . 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 [START_REF] Campi | Wait-and-judge scenario optimization[END_REF].

Theorem 2. [25, Theorem 2] Consider a convex scenario program defined as in [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF]. Under assumptions 1, 2, stationary P and IID samples in D N the solution θ of (16) satisfies

P N [P f (θ ) ≤ (s N )] ≥ 1 -β, (24)
where the reliability (k) = 1t(k) is the unique solution in [0,1] of the following polynomial equation in t for any k = 0, ..., n θ :

β N + 1 N j=k B j (t; k) - N k t N-k = 0. ( 25
)
Here B j (t; k) = j k t j-k is a binomial expansion.

Theorem 2 gives a generalization bound P f (θ ) ≤ (s N ) which holds up to a confidence level 1 -β. The value of the reliability parameter (s N ) is determined from (25) after observing s N = k support constraints which is always ≤ n θ for convex problems. In contrast to the bound in [START_REF] Calafiore | The scenario approach to robust control design[END_REF], (s N ) is a-posteriori computed after enumeration of the set of support scenarios and its cardinality s N in correspondence of θ . An important consequence of this theorem is that a-posteriori observing a solution that has complexity of s N form a dimension n θ > s N does not give equal robustness guarantees as obtaining a solution from a problem in dimension n θ = 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 [START_REF] Campi | Wait-and-judge scenario optimization[END_REF] 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.

Bounds for convex problems with discarded samples

An extension of Theorem 1 allows for k samples δ to be intentionally discarded from D N , for instance, the ones making program [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF] unfeasible. As a result, the optimized θ enjoys an improved J(θ ) for the cost of a weaker certificate of generalization [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF].

Theorem 3. [87, Theorem 2.1] Consider a convex scenario program defined as in Eq. [START_REF] Campi | Interval predictor models: Identification and reliability[END_REF] where k = N 0 scenarios have been intentionally removed from the data set D N . Under assumptions 1, 2 and stationary P and IID samples in D N-k the distribution of P f (θ ) is bounded by:

P m [P f (θ ) ≤ ] ≥ 1 - k + n θ -1 k k+n d -1 i=0 m i i (1 -) m-i = 1 -β, (26)
where n θ is the number of optimization variables and β ∈ [0, 1] is a confidence parameter.

Theorem 3 is used in practice as follows: a desired confidence parameter β and N o 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 θ a reliability parameter (N, β, n θ , N o ) is computed from the right hand side of [START_REF] Ramponi | Expected shortfall: Heuristics and certificates[END_REF]. This gives an upper bound P f (θ ) ≤ (N, β, n θ , N o ) on the violation probability which is a-priori computed, i.e., is obtained before observing θ solution of program 2.

Bounds for soft-constrained convex problems

To extend the scope of scenario-based generalization bounds to the soft-constrained case [START_REF] Vapnik | Statistical Learning Theory[END_REF], 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. [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Given assumptions 1, 2, 3, stationary P and IID samples in D N the probability P f (θ ) is bounded by:

P N [ (k) ≤ P f (θ ) ≤ (k)] ≥ 1 -β, (27)
where k = s N is the number of support constraints and (k) = max{0, 1t(k)}, (k) = 1t(k) and {t(k), t(k)} are solutions of a polynomial equation in t:

B N (t; k) = β 2N N-1 j=k B j (t; k) + β 6N 4N j=N+1 B j (t; k), (28)
where k = 1, ..., N -1 is the number of support constraints. For the special case k = N, the upper bound is set to (k) = 1 and the lower bound is obtained solving

1 = β 6N 4N j=N+1 B j (t; k). (29)
Theorem 4 gives a-posteriori upper and lower bounds on the violation probability of the solution of [START_REF] Vapnik | Statistical Learning Theory[END_REF]. The bounds can be computed in practice from the number of constraints N by selecting a confidence β and enumerating the support constraints s N i.e., [ (N, β, s N ), (N, β, s N )]. Note that, for soft-constrained convex programs, s N = N i=1 1{g(δ i ; θ ) ≥ 0}. The number of support constraints coincides with the number of scenarios for which g(δ i ; θ ) ≥ 0 at the optimum θ .

Bounds for minimax problems

Let us define the ordered statistic of the constraint function values as γ 1 ≥ γ 2 ≥ ... ≥ γ N where the quantity γ i = g(δ i ; θ ) is the i th higher value of the constraint g in correspondence of the optimum θ . The probability of error for an empirical level γ i is given by P f,i = P[δ ∈ ∆ : g(δ; θ ) > γ i ], and the corresponding reliability is R i = P[δ ∈ ∆ : g(δ; θ ) ≤ γ i ]. Note that γ i ≥ γ j for any i < j and, thus, the empirical probabilities of error are P f,i ≤ P f, j for i < j. Theorem 5 prescribes a formal bound on the joint distribution of the probabilities of error P f,n θ +1 , ..., 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 θ +1 , .., P f,N is bounded by

P N [P f,i ≤ n θ +1 , ..., P f,N ≤ N )] ≤ F θ,N ( n θ +1 , ..., N ) = β, (30)
where

F n θ ,N = N! n θ ! n θ +1 0 τ n θ n θ +1 • • • N 0 1 {0≤τ n θ +1 ≤...τ N ≤1} dτ N ...dτ n θ +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 θ +1 ≥ P f,i for i = 1, .., n θ , the bound n θ +1 on the level P f,n θ +1 automatically applies also to P f,i i = 1, .., n θ .

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 γ. This is a non-parametric probability box which bounds the reliability at the empirical violation levels γ. For an overview on probability boxes the interested reviewer is reminded to [START_REF] Ferson | Constructing probability boxes and dempster-shafer structures[END_REF][START_REF] Schöbi | Global sensitivity analysis in the context of imprecise probabilities (p-boxes) using sparse polynomial chaos expansions[END_REF][START_REF] Zhu | A probability-box-based method for propagation of multiple types of epistemic uncertainties and its application on composite structural-acoustic system[END_REF][START_REF] Rocchetta | Assessment of power grid vulnerabilities accounting for stochastic loads and model imprecision[END_REF][START_REF] Rocchetta | A post-contingency power flow emulator for generalized probabilistic risks assessment of power grids[END_REF] and propagation strategies on different applications. The upper bound of the probability box is constructed as follows: [START_REF] Carè | Scenario min-max optimization and the risk of empirical costs[END_REF]:

F(γ) =              1 if γ ≥ γ n θ +1 1 -i if γ i+1 < γ ≤ γ i , (i = n θ+1 , ..., N -1) 1 -N if γ < γ N , (31)
whilst the lower bound is given by,

F(γ) =              1 -1 if γ ≥ γ 1 , 1 -i if γ i ≤ γ ≤ γ i-1 , (i = 2, ..., N) 0 if γ < γ N , (32)
where i and i are obtained for a confidence level β from Theorem 4. The function F(γ) and F(γ) are the upper and lower CDFs bounding empirical reliability levels

F(γ) ≤ P[δ ∈ ∆ : g(δ; θ ) ≤ γ] ≤ F(γ),
with a confidence at least 1 -2β.

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 (θ ) is bounded by,

(N, β, s N ) =            1 if s N = N 1 -β N( N s N ) 1 N-s N otherwise , (33)
where s N is the number of support scenarios and β ∈ [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 [START_REF] Calafiore | A learning theory approach to the construction of predictor models[END_REF], i.e., programs for which g(δ; θ) or J(θ) are non-convex functions of θ. 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 the remainder of this paper we will focus on convex programs and reliability bounds for convex problems.

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 β 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 θ 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 θ before solving the optimization and observing the optimized solution, i.e., a-priori [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF][START_REF] Campi | A sampling-and-discarding approach to chance-constrained optimization: Feasibility and optimality[END_REF]. 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 θ .

In contrast, a-posteriori evaluation of the solution's complexity can provide a tighter quantification of the epistemic uncertainty affecting its reliability [START_REF] Campi | Wait-and-judge scenario optimization[END_REF][START_REF] Garatti | Complexity-based modulation of the data-set in scenario optimization[END_REF][START_REF] Garatti | Risk and complexity in scenario optimization[END_REF]. 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 θ , the number s N is only available after optimizing θ and is a random variable that depends on the multi-samples extraction 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 () 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 θ (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 β = 10 -6 , assuming n θ = 100, s N = 95 and N o = 5 discarded samples. For additional details and a numerical comparison of the bounds see, for instance, [START_REF] Geng | Data-driven decision making with probabilistic guarantees (part 1): A schematic overview of chance-constrained optimization[END_REF].

(N, β, n θ ) 1 [24] Hard-constrained 1 (N, β, s N ) 2 (N, β, s N ) 6 [25] [28] Hard-constrained (discarded) 2 (N, β, n θ , N o ) 3 [87] Hard-constrained (discarded) 2 (N, β, s N + N o ) 2 [25] Soft-constrained 3 [] Soft-constrained 3 (N, β, s N ), (N, β, s N ) 4 [54] Minimax 4 (N, β, n θ ) 1 [24] [94] Minimax 4 1 (N, β, n θ ), ..., N (N, β, n θ ) 5 (N, β, s N ) 6 [31] [28]
Table 2 presents the bounds as functions = (N, β, ...), 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 [START_REF] Campi | Interval predictor models: Identification and reliability[END_REF]:

• Design of an identification experiment, N = N( , β, n θ , N o ). The problem of finding the minimum number of observations N which are necessary to obtain a certain reliability score given β, n θ and N o .

• Maximal complexity and model selection, n θ = n θ (N, , β, N o ). The problem of determining the maximum number of the fitting coefficients n θ to achieve a reliability at least given N samples, N o outliers, and a confidence β. 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 θ ≤ 7 design variables, we can ensure a-priori a reliability level of at least 0.99.

• Data cleansing, N o = N o (N, , β, n θ ). 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 ?"

Figure 1: A comparison of the scenario bounds summarized in Table 2.

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(δ; θ) ≤ 0 on the lower and upper bounds. A relaxation of the constraints enables the condition g(δ; θ) ≤ 0 to be violated, i.e., some of the data points can fall outside from the optimized bounds if this yield a significant improvement of the area metric. Consider the following soft-constrained IPM training program.

IPM program 4. (Soft-Constrained IPM design) θ , ζ = arg min θ, ζ i ≥0 A(θ) + ρ N i=1 ζ i : s.t. f l (x i ; θ) -y i ≤ ζ i , i = 1, ..., N, y i -f u (x i ; θ) ≤ ζ i , i = 1, ..., N, ( 34
)
where ζ is a N-dimensional vector of non-negative slack variables relaxing the samples constraints on the upper and lower bounds.

If ζ i > 0 at the optimum, the sample (x i , y i ) will fall outside the bounds defined by I N (x; θ ). Conversely, if ζ 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 ∈ I N (x i ; θ ). Note that if the cost parameter ρ → ∞ 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 A(θ ), 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 ζ leading to a decision space of dimension N + n θ > N which can therefore lead to a solution having complexity s N > n θ . 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:

P N ≤ P[(x, y) ∈ ∆ : y I N (x; θ )] ≤ ≥ 1 -β.
By selecting a confidence level, an epistemic interval on the probability of error is derived as follows:

P f (θ ) ∈ [ (N, β, s N ), (N, β, s N )].
Similarly, an interval on the IPM reliability holds with the same confidence level and is given by:

R(θ ) ∈ [1 -(N, β, s N ), 1 -(N, β, s N )].
The program in Eq. ( 34) is convex in θ and, thus, s N can be easily obtained from the number of violating scenarios (for which ζ 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 A. For convex programs, a support constraint is always active at the optimum2 and a conservative estimate of s N is given by s N = |S| where the set of support constraints is

S = (x i , y i ) ∈ D N : ζ i > 0, ∨ y i -f u (x i ; θ ) = 0, ∨ f l (x i ; θ ) -y i = 0 ,
where ∨ is the logic disjunction defining the 'or' boolean operator.

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 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.

Case study I

The DGM of the first case study is defined as follows [START_REF] Crespo | Interval predictor models with a formal characterization of uncertainty and reliability[END_REF][START_REF] Lacerda | Interval predictor models for data with measurement uncertainty[END_REF]:

y(x; σ) = x 2 cos(x) -e -x 2 sin(3x) -x -cos(x 2 ) + xN(0, σ) (35)
where x ∈ X is an input variable in X = [-10, 5], and N(0, σ) is a Gaussian distributed random noise with zero mean and standard deviation σ. A data set D N is obtained by uniform sampling of N realization of x i in X and computing y i with a randomized noise value.

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

A data set D N=150 is sampled from Eq. ( 35) with σ = 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 ψ(x) and n l = n u fitting coefficients ranging between 3 and 20. Thus, the number of optimization variables n θ 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 θ 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 θ . Monte Carlo estimators of the true probability of error 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 θ 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 θ = 11 presented in Figure 2, which render A = 0.79 and 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 3). The a-posteriori and a-priori Scenario bounds from Theorems 1 and 2, are also presented for a confidence parameter β = 10 -6 (almost certainty) and verified via Monte Carlo (MC) estimation of the true failure probability.

n l = n u 3
(10, 150, β) = 0.23 for n l = 5 and (16, 150, β) = 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 θ the a-posteriori complexity is generally s 150 < n θ , and sometimes much lower than n θ . 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.

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 [START_REF] Shafer | A mathematical theory of evidence turns 40[END_REF] 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 -= 0.72 (see Table 3) whilst the IPM obtained discarding samples is only R(θ ) > 1 -= 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.

Minimax layer

The Minimax layer program ( 10) is solved using a data set of N = 2000 samples and a regression function f (x; θ) defined by radial basis functions with scale parameter ν = 1 and n θ = 20 fitting coefficients. Figure 4 presents the results. The optimized f (x; θ ) is displayed by the red curve on the left panel and it has been used to construct two IPMss [ f (x; θ ) -γ 1 , f (x; θ ) + γ 1 ], and [ f (x; θ ) -γ 10 , f (x; θ ) + γ 10 ], for the empirical levels γ 1 and γ 10 , respectively. The empirical level γ 1 = max i=1,..,N |y if (x; θ )|, 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 γ 10 corresponds to the 10 th highest difference |y if (x; θ )|. 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 γ 1 (solid lines in the left panel) and γ 10 (dotted lines in the left panel). The accuracy of the predictor [ f (x; θ ) ± γ 1 ] is given by an area 2γ 1 which results in a larger area when compared to the area of the second predictor 2γ 10 . Theorem 5 guarantees that the reliability of the predictor [ f (x; θ ) ± γ 1 ] is at least R(θ ) > 0.99. Although much more accurate, the reliability for a level γ 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 [START_REF] Dubois | Possibility theory, probability theory and multiple-valued logics: A clarification[END_REF] 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; θ ) 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.

Results for the proposed soft-constrained program

Figure 5 presents two examples of IPMs trained using program [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF]. Radial basis functions ψ(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 ρ = 0.5 is presented whilst the bottom panel presents an IPM trained with softer-constraints ρ = 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 ρ worsen the accuracy of the IPM leading predictive intervals which are (in average) wider. Specifically, a ρ = 0.5 leads to a A = 0.582 whilst the area for a ρ = 0.05 is only A = 0.439.

The IPM trained with hard constraints is more reliable. Specifically, we can guaranteed that the reliability of the IPM trained for a ρ = 0.5 and ρ = 0.05 lay, respectively, in R ∈ [0.96, 1.0] and R ∈ [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 R = 0.992 and R = 0.956 for ρ = 0.5 and ρ = 0.05, respectively. For both cases, R(θ ) ∈ [1 -, 1 -] as prescribed by Scenario theory.

The coefficients of soft-constrained radial basis bounds of degree n θ = 20 are optimized using program [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF] and for 50 values of ρ between 0.01 and 0.5. We examine 20 training sets D N of size N ∈ [50, 10 3 ]. Figures 6 and7 present the results of this analysis. In Figure 6 the area A, the scenario bounds [ , ] 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 ρ value). We can observe a trade off between area and both [ , ] and P f . Also, the failure probaiblity is positively correlated to and 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 P f --. 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 (blue markers). In contrast, smaller N leads to a generalization interval which is tighter from above, i.e., P f results closer to (larger red markers).

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 noncollocated sensor/actuator presented in the NASA langley challenge on optimization under uncertainty [START_REF] Crespo | The nasa langley challenge on optimization under uncertainty[END_REF] on y-axis against the 20 values of N (x-axis) and the 50 values of ρ (iso-N markers ). Interestingly, a higher N leads to a P f that is closer to the lower bound 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 [START_REF] Crespo | The nasa langley challenge on optimization under uncertainty[END_REF]. 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 [START_REF] Crespo | The nasa langley challenge on optimization under uncertainty[END_REF]. 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)} T t=1 100 i=1 and the total number of sample pairs is 5 × 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 [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF] to define 10 4 constraints and identify an optimal IPM. A cost parameter ρ = 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 ν = 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 R = 0.9906 and support the validity of the epistemic bounds.

Note that R is obtained averaging over the N = 5 × 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 P[z 2 (t) ∈ I N (t; θ ), ∀, t = 1, ..., T ]. that is the probability of I N (t; θ ) 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.

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 [START_REF] Rocchetta | On-line bayesian model updating for structural health monitoring[END_REF] for more details on the model). For reproducibility sake, we define the DGM of this The results for the Minimax layer IPM defined by the 20 st highest empirical level γ 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 minutes3 . Figure 10 presents the bounds I N (x; θ ) optimized via program (3), program (8) optimally discarding N o = 400 samples, and the soft-constrained program [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF] with ρ = 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.

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. Scenariobased bounds hold distribution-free, non-asymptotically and their width depends on sample size available to characterize the uncertain process. Moreover, we show that simpler models generally lead to 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 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. [START_REF] Faes | On the robust estimation of small failure probabilities for strong nonlinear models[END_REF] for preliminary comparison between Kirging 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 [START_REF] Nesterov | Subgradient methods for huge-scale optimization problems[END_REF][START_REF] Lancia | Large-Scale Linear Programming[END_REF].

Investigate non-convex IPM optimization programs and new epistemic bounds for non-convex cases. In particular, 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 formally 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.

Definition 2 .

 2 (Non-reducible support set) A support set S ⊆ D N is a k-tuple S = {δ i 1 , ..., δ i k } for which the solutions of a scenario program with constraints defined by D N is identical to the same program with S in place of D N . A set S ⊆ S is non-reducible if, for any δ ∈ S , the constraints defined by S \ δ 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 | • | 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 ∈ N 0 and scenario set D N , the optimal solution θ of a scenario program exists and is unique.

Figure 2 :

 2 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.

Figure 3 :

 3 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.

Figure 4 :

 4 Figure 4: In the panel on the left: the optimized regression function f (x; θ ) and two minimax layer IPMs constructed for the empirical levels γ 1 and γ 10 . In the panel on the right: the probability box bounding the reliability of minimax layer IPMs for any γ i , i = 1, ..., N. The P-box holds with a confidence β = 1 -2 • 10 -4 .

Figure 5 :

 5 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 σ = 3. The red x markers show the failed validation samples. The top and bottom panels display the optimized IPM resulting from cost parameters ρ = 0.5 and ρ = 0.05, respectively.

Figure 7 :

 7 Figure

Figure 10 : 4 Table 5 :

 1045 Figure 10: A comparison of radial basis IPMs trained using programs (3), discarding N o = 400 samples, and the soft-constrained program with ρ = 0.01.Program Eq.[min]A(θ ) R(θ ) [1 -, 1 -] Theor. (3) 3.4 0.72 0.924 [0.881, 1.0] 2 (8) N o = 200 23.4 0.59 0.859 [0.724, 1.0] 3 (8) N o = 400 19.3 0.44 0.731 [0.519, 1.0] 3 (10) γ 20 0.1 3.30 0.989 [0.973, 1.0] 5 (34) ρ = 0.05 12.0 0.62 0.874 [0.824, 0.91] 4 (34) ρ = 0.01 22.5 0.52 0.786 [0.726, 0.83] 4Table 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.

Table 2 :

 2 Summary of the reviewed reliability bounds 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 () 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.

		Scenario	convex	non-convex	Type
	Type	Program	Bounds	Theorem Bounds	Theorem a-prio. a-post.	Ref.
	Hard-constrained	1

Table 3 :

 3 Area metrics, computational time and number of support constraints for N = 150 and IPM trained with program (

			4	5	7	9	11	15	20
	CPU Time [s]	0.02	0.1	0.3	1.7	6.3	17.1	49.5	111.7
				Accuracy
	A(θ )	0.85	0.84	0.79	0.66	0.54	0.42	0.33	0.32
			Epistemic bounds on P f
	n θ	6	8	10	14	18	22	30	40
	(n θ , N, β) a-priori	0.158 0.18	0.20	0.24	0.28	0.31	0.37	0.45
	6 (s N , N, β) a-posteriori 0.19 s N	8 0.21	10 0.23	14 0.27	15 0.28	16 0.29	19 0.32	19 0.32
				Validation
	P f (θ ), MC	0.06	0.072 0.065 0.105 0.151 0.153 0.215 0.230

Table 4 :

 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 β = 10 -6 .

	N o	0	10	50	100	125	150	175	200	225
	CPU Time [s]	1.19	1.090 0.87	0.90	0.91	98.40 90.99 1.05	85.08
				Accuracy
	A(θ )	0.56	0.52	0.45	0.43	0.40	0.38	0.39	0.39	0.40
			Epistemic bounds on P f
	(n θ , N, β, N o ) a-prio.	0.04	0.07	0.12	0.172 0.195 0.21	0.24	-	-
	20 (s N + N o , N, β) a-post. 0.028 s N	32 0.037 0.065 0.096 0.110 0.128 0.139 0.156 0.17 70 122 146 176 196 224 251
				Validation
	P f (θ ) MC	0.0127 0.018 0.041 0.060 0.077 0.092 0.097 0.114 0.123

  . The 18

	softer constraints: larger failure probability
	but narrow predictive interval
	harder constraints: small failure probability
	but wider predictive intervals
	softer constraints
	harder constraints

compromise Figure 6: The results for 1000 IPMs trained with program (34). The 3 marker types indicate different levels of relaxation for the constraints.

1+exp(-c j •x i )

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).

Note that an active constraint is not always a support constraints 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.

Note that the computational cost can be reduced of an order of magnitude by adopting linear programming solvers

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)

process by estimating the means and standard deviations of the frequency response over a set knots on the support X. Figure 9 presents the mean values (black marked line) and mean ± the standard deviations (blue dashed and dotted lines). A data set D N is generated by sampling x i ∈ 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 of 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.

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 θ achieve better results. Thus, we define the bounds by linear combinations of n l = n u = 80 radial basis and fitting coefficients θ. A scale parameter ν = 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 A(θ ) = 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(θ ) > 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 θ , can be selected (however this likely leads to a lower accuracy).

If the average width of the IPM is too large, 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 [START_REF] Shafer | A mathematical theory of evidence turns 40[END_REF] 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 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(θ ) > 1 -(N, β, s N + 200) = 0.724 when 200 are discarded, and R(θ ) > 1 -(N, β, 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 [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF]. For a costs of violation ρ = 0.05 and ρ = 0.01, [START_REF] Carè | The wait-and-judge scenario approach applied to antenna array design[END_REF] 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 ρ increases the bounds on the reliability get worse. This is due to the lower cost of constraint violations leading to larger s N .