Conradus Van Mierlo

Augustin Persoons

Matthias G R Faes

David Moens

Robust design optimisation under lack-of-knowledge uncertainty

Keywords: robust optimisation, interval analysis, Gaussian Process modeling, efficient global optimisation

Design optimization is common practice in engineering where the goal is to find the optimal combination of design parameters under prescribed constraints. However, some parameters may be impossible to define in a deterministic sense and may only be known with significant uncertainty. This limitation has led to an alternative definition of design optimality called robustness, where attention is payed to the variation around the optimal performance. Straightforward methods to solve robust optimization problems are usually limited in two ways: (1) the computation burden of the so-called 'double-loop' optimization problem hinders application to realistic models, and (2) the formalisms are typically limited to probabilistic descriptions of the uncertainty. This paper presents a formulation of the robust optimization problem under interval uncertainty and proposes a new approach taking advantage of the so-called adaptive Gaussian processes to solve it efficiently. The proposed surrogate approach mitigates the computational burden of the resolution, and a dedicated learning function is proposed to ensure iterative minimization of the surrogate modelling error and convergence towards the robust optimum. The algorithm uses a stopping criterion related to the level of confidence associated with the optimality of the solution. The approach is illustrated on six analytical and engineering benchmark problems.

Introduction

Current engineering practice involves the development and design of products that span an ever growing field of applications, while the performance of these products should also be guaranteed under a wide range of circumstances. In other words, the performance of a product should be only minimally affected by, e.g., load variations, changing environments, boundary conditions.

The idea of products and processes that are insensitive to variations, e.g., in manufacturing, was pioneered by Genichi Taguchi who first applied his methodology on electrical circuits [START_REF] Taguchi | Quality engineering (taguchi methods) for the development of electronic circuit technology[END_REF][START_REF] Taguchi | Performance analysis design[END_REF]. However, the description of these variations, including the details about their underlying probability density functions (PDF's), is in general a challenging task. The main reasons for this are that the corresponding quantities are inherently variable, e.g. wind loads, there is incomplete knowledge about the quantity, e.g. direct measurement is challenging, or the designer is faced with a combination of both [START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF]. Additionally, in an early design stage, where the fundamental design decisions are made, only rough estimations of the quantities influencing the performance might exists. Historically, in engineering practice uncertainties are covered by safety factors. Although this approach is very straightforward, these safety factors will not provide information about the actual conservatism in the design. Therefore, numerous techniques for uncertainty quantification have been introduced during the last decades to account for these uncertainties. Typically, these techniques are categorised as probabilistic and possibilistic approaches [START_REF] Stefanou | The stochastic finite element method: Past, present and future[END_REF]. The latter includes techniques as: interval [START_REF] Faes | A multivariate interval approach for inverse uncertainty quantification with limited experimental data[END_REF][START_REF] Moore | Interval analysis[END_REF], fuzzy sets [START_REF] Hanss | Applied fuzzy arithmetic[END_REF], information gap methods [START_REF] Ben-Haim | Info-gap decision theory: decisions under severe uncertainty[END_REF], and imprecise probabilities [START_REF] Faes | Engineering analysis with probability boxes: a review on computational methods[END_REF][START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF]. In general, probabilistic methods are best suited for aleatory uncertainties as they describe non-determinism via random variables defined by their joint probability distributions, while possibilistic approaches are usually well suited to cover both aleatory and epistemic uncertainties.

In addition to the variety of possibilistic methods, different definitions of the robustness are proposed in literature; the relevance of which depends on, e.g., the application and the available information. For a review of different robustness measures under probabilistic uncertainty the reader is referred to the work in [START_REF] Moritz Göhler | Robustness Metrics: Consolidating the Multiple Approaches to Quantify Robustness[END_REF][START_REF] Zhang | Nuanced robustness analysis with limited information[END_REF][START_REF] Park | Robust design: an overview[END_REF]. In the context of possibilistic uncertainties, robustness definitions have been introduced in the framework of information gap theory [START_REF] Kuczkowiak | Experimental Validation of an Info-Gap Uncertainty Model for a Robustness Analysis of Structural Responses[END_REF][START_REF] Roach | Comparison of robust optimization and info-gap methods for water resource management under deep uncertainty[END_REF], convex models [START_REF] Au | Robust design of structures using convex models[END_REF], and for fuzzy sets [START_REF] Beer | Designing robust structures-a nonlinear simulation based approach[END_REF]. The definitions in these works are mainly based on two criteria: the first is minimising the variation of the output [START_REF] Beer | Designing robust structures-a nonlinear simulation based approach[END_REF], and the second is to optimise simultaneously both the output (e.g. performance) and its variance around the optimal value [START_REF] Zang | A review of robust optimal design and its application in dynamics[END_REF][START_REF] Hu | An efficient robust optimization method with random and interval uncertainties[END_REF][START_REF] Mitseas | Robust design optimization of structural systems under evolutionary stochastic seismic excitation[END_REF]. In addition to the definition of robustness, a range of methods have been developed for its evaluation, with sampling strategies for most mixed uncertainty problems [START_REF] Beer | Designing robust structures-a nonlinear simulation based approach[END_REF][START_REF] Lee | Robust design optimisation using multiobjective evolutionary algorithms[END_REF], forward or inverse propagation [START_REF] Taguchi | Quality engineering (taguchi methods) for the development of electronic circuit technology[END_REF], meta-model assisted methods [START_REF] Ur Rehman | Efficient kriging-based robust optimization of unconstrained problems[END_REF], and fully decoupled methods for reliability based design optimisation [START_REF] Faes | Fully decoupled reliability-based design optimization of structural systems subject to uncertain loads[END_REF].

This work focuses on developing a meta-model assisted method to determine the robustness at different design points. The meta-model that is used is a Gaussian Process (GP) model also known as Kriging [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF], which is used in this contex as an emulator of the physical model. After calibration of the GP-model on a set of evaluated points, i.e., Design of Experiments (DOE), the model is fast to evaluate. Based on this easy to evaluate GP model, fast approximations can be made about the underlying problem, i.e., numerical model, and this approximation can be improved by increasing the calibration points in the DOE. The well-known framework of Efficient Global Optimisation (EGO) [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] successfully exploits the GP mean and variance to select additional calibration points and improve on the predicted minimum. In the specific case of interval uncertainties, the GP is used to estimate the interval width in un-sampled regions, including the confidence bounds about this estimate. Hence, the GP estimate can be used in place of the actual model for the optimization problem. The estimation will be affected by a modeling error but can be bounded by a confidence interval. An improvement function is proposed that finds the next point to evaluate as a compromise between its estimated robustness and the uncertainty regarding its estimation (high GP variance). The improvement function in this paper is based on the work of M. De Munck et al. [START_REF] De Munck | An efficient response surface based optimisation method for non-deterministic harmonic and transient dynamic analysis[END_REF]. However, in this work some adaptions are proposed to the improvement function to efficiently perform the robust optimisation. To solve the robustness optimisation efficiently two improvement functions are combined: first an improvement of the interval width throughout the domain and second an improvement towards the most robust design. The combination of these two improvement functions provides a powerful improvement function that refines the GP model both globally and locally around the most robust design point.

The proposed Robustness under Lack-of-Knowledge method is abbreviated as RuLoK. This paper is structured as follows: Section 2 describes the terminology and notation that is used. In section 3 the robustness measure under interval uncertainty is introduced, while Section 4 provides the details towards the Gaussian Process model that is used. Section 5 describes the adaptive sampling strategy that is used to calibrate the Gaussian process and the performance of this method is demonstrated on a number of cases in Section 6. Finally, in section 7 a discussion about the results is held before conclusions are drawn in Section 8.

Terminology and notation

In this paper, a vector is indicated as lower-case boldface characters x, matrices are expressed as upper-case boldface characters X and interval parameters are indicated using apex I: x I . Further, a distinction is made between design parameters and uncertain parameters.

Design parameters: z ∈ Z ⊆ R nz with Z the set of admissible designs and n z ∈ N, are controlled and part of the design problem, e.g., plate thickness, hole diameters.

Uncertain parameters: x ∈ x I ⊆ IR nx are uncontrollable, purely epistemic parameters with n z ∈ N, typically modeled by intervals that represent, e.g., wind loads, electrical resistance, transmission parameters. The uncertain parameter vector x is represented as an interval

vector x I = [x I 1 , x I 2 , . . . , x I nx ]
, with x I i , i, . . . , n x the i th parameter interval. An interval is considered closed when both the upper and lower bounds are a member of the interval. The domain of a real-valued interval is denoted as IR.

Robustness under lack-of-knowledge uncertainty

The uncertainty considered in this work is purely epistemic in nature and results from a lackof-knowledge about the exact value of the parameter. In practice, this kind of uncertainty is encountered when the best estimate of a parameter is limited to a range of possible values, even when its based on all available data and/or knowledge. The real value of the quantity, be it deterministic or variable, is in this case represented by the bounds between which it is deemed to lie. Precisely, an interval is defined as:

x I = [x; x] = {x ∈ R nx | x ≤ x ≤ x}, (1)
where x denotes the lower bound and x denotes the upper bound. In addition, an interval can be represented by the centre point x = x+x 2 and radius ∆x = x-x 2 of the interval.

Propagation of interval valued uncertainty

In this work the model m is a continuous function on R, which is parameterised by a parameter vector θ θ θ. The parameter vector consists out of two parts θ θ θ = {x, z}, with x the uncertain parameters and z the design parameters. The number of elements in the parameter vector are indicated by n θ θ θ = n x +n z . By solving the model m the parameter vector θ θ θ is transformed R n θ θ θ → R to a scalar response quantity y ∈ Y ⊂ R, with the set of admissible model parameters Y, defined as:

m : y = m(θ θ θ). (2)
The main goal of the interval analysis is to identify the extremes of the set of system responses ỹ. Since finding the set ỹ is in general computationally intractable, the exact solution set is often approximated by a realisation set ỹs defined as [START_REF] Faes | Recent trends in the modeling and quantification of non-probabilistic uncertainty[END_REF]:

ỹs = y j | y j = m(θ θ θ j ); x j ∈ x I ; j = 1, . . . , n q . ( 3
)
The set ỹs is typically constructed by performing n q deterministic evaluations y j = m(θ θ θ j ) of the numerical model, with y j the response of the j th solution. For each of these n q solutions, a sample is taken within the range of the interval x I . The main challenge herein is choosing x j such that ỹs is a accurate approximation of ỹ. A first way to obtain such approximation is to follow an optimisation approach. Here, the exact solution set ỹ is approximated by an accurate interval for the one dimensional case. For the higher dimensional case a conservative approximation is made about the hyper-cubic solution set in higher dimensions y I = [y I 1 , y I 2 , . . . , y I ny ], with ỹ ⊆ y I . The corresponding optimisation problem is defined as:

y = min x∈x I m(θ θ θ), y = max x∈x I m(θ θ θ), (4)
where y I = [y; y] is the solution interval. When a global minimum or maximum is found through optimisation, the exact output set bounds are obtained. However, it should be noted that the behaviour of the goal function with respect to the uncertain parameters is unpredictable in the case of strongly non-linear problems, which makes the computational effort highly problem dependent [START_REF] Moens | Non-probabilistic finite element analysis for parametric uncertainty treatment in applied mechanics: Recent advances[END_REF].

There is a special case for monotonic problems, where the vertices of the hyper-cubic input space are sampled, called the vertex method, introduced by Dong and Shah [START_REF] Dong | Vertex method for computing functions of fuzzy variables[END_REF]. Following this method the output set is determined exactly within 2 nx evaluations. However, the underlying assumption is that the model output behaves monotonically with respect to the input parameters, which is not true in general. Other approaches are intrusive methods to solve interval problems, which have been proposed in [START_REF] Muscolino | Stochastic analysis of structures with uncertain-but-bounded parameters via improved interval analysis[END_REF], and interval arithmetic methods as proposed in [START_REF] Sofi | A novel interval finite element method based on the improved interval analysis[END_REF].

Defining robustness in the case of interval valued uncertainty

As mentioned in the introduction, multiple definitions of robustness exist, depending on the context and application. In this work, it is proposed to define robustness as the design with minimum variation in the performance given a well-defined input uncertainty. Following this definition, robustness can be defined as the ratio of input uncertainty to the output uncertainty. However, quantifying this uncertainty is non-trivial in general. Therefore, the focus lies on the interval radius as a measure for the uncertainty. In this way, this robustness measure can be regarded as an interval counterpart to robustness measures that minimize the variance of the performance. For a case with one interval valued input parameter, the input and output uncertainty are represented respectively by the scalar interval radius ∆x and the associated scalar output interval radius ∆y.

The output radius is a function of the design parameter z and should be evaluated for multiple designs z ∈ Z. The robustness for this case is defined as:

R(z) = ∆x ∆y(z) = x -x y(z) -y(z) , (5)
Since ∆x is independent of the design z, finding the most robust design z * is reformulated to the minimisation of the output uncertainty, which can be evaluated for multidimensional cases, defined by:

z * = argmin z∈Z [y -y] = argmin z∈Z [max x∈x I m(θ θ θ) -min x∈x I m(θ θ θ)]. (6)
Figure 1 illustrates the proposed robustness measure R for a point z * and shows the associated upper bound y(z) and lower bound y(z), in red and blue. The point z * is also the point with the maximum robustness R, indicated in orange. As suggested from Equation ( 6) finding the robustness of just one design involves a global optimisation to construct the conservative approximation of the solution set ỹ, which should be repeated for each of the design points in Z. Thus, crude optimisation of the problem described in Equation ( 6) involves two other optimisation problems: first an optimisation that actively looks for the upper-bound y, and second, an optimisation that searches the lower bound y, both for the same design z. Therefore, crude optimisation is a time consuming effort, as this would involve a large number of evaluations of the model m under consideration. In an attempt to alleviate this problem, the next section discusses the use of a well-designed Gaussian process model G that could be used in place of the model m.

y

z i y(z i ) y(z i ) z * i R(z * i )
Figure 1: Illustration of the optimal robust design points R(z * ) (orange) for the upper and lower bounds y and y for a specific design parameter z i .

Gaussian process model for robustness under interval uncertainty

This section provides a short theoretical summary of Gaussian Process (GP) models or Kriging [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF][33], an introduction with examples is also available in [START_REF] Rasmussen | Gaussian processes in machine learning[END_REF]. A GP model is a stochastic meta-model that assumes m(θ θ θ) to be a realisation of a Gaussian process, which is defined as [START_REF] Santner | The design and analysis of computer experiments[END_REF]:

G = β T f (α) + σ 2 (x, Ω), (7)
with the first term being a deterministic regression model with

f (α) = {f 1 (α), . . . , f k (α)} a set
of arbitrary basis functions, and β T a vector of regression coefficients. The second term consists of a zero-mean, unit variance, stationary Gaussian process (x, Ω) scaled with a constant variance of the Gaussian process σ 2 . The underlying probability space of the Gaussian process is represented by Ω and the correlation between two points r and r is defined by the covariance function K(r, r , l c ), with l c the characteristic length or other hyper-parameters. In general, one refers to the covariance matrix K where the covariance is determined for all points in a domain.

The reader may refer to [START_REF] Alvarez | Kernels for vector-valued functions: A review[END_REF] for details about different covariance functions in Gaussian processes.

In this paper two well-known covariance functions are used: The Gaussian kernel (also known as squared-exponential covariance function) and the Matérn 5 2 kernel.

The GP-model is then calibrated on an initial design of experiments x DOE obtained from, i.e., Latin hyper-cube sampling and their observed results y DOE . Conditional on the observed data the mean and the variance of the Gaussian process can be estimated [START_REF] Santner | The design and analysis of computer experiments[END_REF]:

µ gp = f (x) T β + r(x) T K -1 (y DOE -F β), (8)
σ 2 gp = σ 2 1 -r T (x) + u T (x)(F T K -1 F) -1 u(x) , (9)
with F the matrix of the observed trend, r(x) a vector of cross-correlations between predicted points x and observed points, and with:

β = (F T K -1 F) -1 F T K -1 y DOE , (10)
the general least-squares estimate of β and

u(x) = F T K -1 r(x) -f (x). (11)
Equations ( 8) and ( 9) are referred to as the mean and variance of the GP predictor, respectively.

The GP that is used in this work is an interpolating GP, which means that the prediction of the variance at an experimental point x ∈ x DOE tends to zero.

Predicting interval bounds with a Gaussian Process model

In this work a GP-model is used to predict the output of the model m with as input θ θ θ the set of uncertain and design parameters. To this end, µ gp is considered to be the best GP-estimate and σ gp is the confidence over this estimate. For the specific application of estimating an output interval based on the GP-model the main interest goes to the maximum and the minimum response over the complete range of uncertainty. Therefore, the bounds of the response are estimated by:

y gp (z) = µ gp (z) = max x∈x I µ gp (θ θ θ), (12)
y gp (z) = µ gp (z) = min x∈x I µ gp (θ θ θ). (13)
A similar approach can be taken to identify the maximum and minimum of the confidence bounds:

δ µ+σ (z) = max x∈x I (µ gp (θ θ θ) + cσ gp (θ θ θ)), (14)
δ µ+σ (z) = min x∈x I (µ gp (θ θ θ) + cσ gp (θ θ θ)), (15)
δ µ-σ (z) = max x∈x I (µ gp (θ θ θ) -cσ gp (θ θ θ)), (16)
δ µ-σ (z) = min x∈x I (µ gp (θ θ θ) -cσ gp (θ θ θ)), (17)
with cσ confidence bounds. The bounds of the response are estimated for each design point z, based on Equations [START_REF] Zhang | Nuanced robustness analysis with limited information[END_REF][START_REF] Park | Robust design: an overview[END_REF][START_REF] Kuczkowiak | Experimental Validation of an Info-Gap Uncertainty Model for a Robustness Analysis of Structural Responses[END_REF][START_REF] Roach | Comparison of robust optimization and info-gap methods for water resource management under deep uncertainty[END_REF][START_REF] Au | Robust design of structures using convex models[END_REF][START_REF] Beer | Designing robust structures-a nonlinear simulation based approach[END_REF]. Figure 2a illustrates a simplification of the GP-model output for a single uncertain parameter x ∈ x I and a single design variable z ∈ Z. The upper bound of the output interval is determined by Eq. ( 12), indicated by the red line, and using Eq. ( 13) the lower bound is found, indicated by the blue line. In addition, the bounds based on the mean plus variance µ gp + cσ gp are predicted by Equations ( 14) and ( 15), indicated by the red dotted and dashed lines. Similarly, the bounds based on the mean minus the variance µ gp -cσ gp are given by Equations ( 16) and ( 17) are indicated by the blue dotted and dashed lines. Moreover, two designs z * gp and z pot gp are shown, illustrating the predicted behaviour along the uncertain parameter x. Note that in general, for one specific design, e.g., z * gp , the location of x for the predicted upper bound y gp (z * gp ) and the location of x of the maximum of the CI for the upper bound δ µ+σ are different.

z p o t g p = 2 z * g p = 4 1 2 3 4 5 1 2 3 4 5 z x y µ gp µ gp + cσ gp µ gp -cσ gp y gp (z) y gp (z) δ µ-σ (z) δ µ-σ (z) δ µ+σ (z) δ µ+σ (z) (
y gp (z) δ µ-σ (z) δ µ-σ (z) δ µ+σ (z) δ µ+σ (z) y gp (z) ± cσ gp (z) y gp (z) ± cσ gp (z) (
In the second illustration, Figure 2b, the estimated interval bounds are shown by the red and blue lines, with the CI about these estimates indicated by the red area for the upper bound, and blue area for the lower bound. Note that the red area is drawn between the upper bound of the minimum prediction and the upper bound of the maximum prediction by the GP-model. In addition, two designs z * gp = 2 and z pot gp = 4 are highlighted to illustrate the robustness measure.

The robustness in Eq. ( 5) can be calculated based on these bounds given by the GP-model.

Specifically, for the design z * gp the robustness is given by:

R(z * gp ) = x -x y(z * gp ) -y(z * gp ) (18)
with y(z * gp ) -y(z * gp ) the estimated interval width, which corresponds to 2∆y(z * gp ). Moreover, based on the CI it is also possible to estimate the potential interval width for z pot gp , which would potentially have a higher robustness. To make this estimate the confidence bounds about the mean prediction are used:

R(z pot gp ) = x -x δ µ-σ (z pot gp ) -δ µ+σ (z pot gp ) (19)
with δ µ-σ (z pot gp ) -δ µ+σ (z pot gp ) the estimated interval width, which corresponds to 2∆δ gp (z pot gp ). The difference between these two robustness measures is that R(z * gp ) is estimated on the mean and the potential robustness R(z pot gp ) is estimated using the CI. Hence, the learning function introduced in Section 5 will exploit this difference, to search for designs with a potential higher robustness.

Note that changing the constant c in Equations ( 14) to [START_REF] Beer | Designing robust structures-a nonlinear simulation based approach[END_REF] from, e.g., 2σ to 3σ will enlarge the distance between red and blue surfaces.

In general, identifying the minimum and maximum as stated in Equations ( 12) until ( 17) is not trivial and involves numerous calls to the GP model. In addition, the mean and variance of the GP model are hard to use for optimisation as in the general case the problem is non-convex.

Therefore, using a GP model is challenging for global optimisation methods. However, a number of successful strategies have been proposed to efficiently optimise such problems e.g., using branch and bound algorithms as proposed in [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. In this work, the continuous problem is discretised over a grid with a fixed number of points. In that case, the complex problem of identifying the maximum and minimum reduces to identifying the highest value in a set of candidates. Note that this only works efficiently with a low number of parameters, as the computational burden increases exponentially O(n d ) with the d-dimensions of the problem for a full grid. In addition, an associated disadvantage is the finite accuracy achievable by the discretisation of the problem, with a finer discretisation causing a higher computational burden. The effects of discretisation can be mitigated in low dimensional problems by using a high number of grid points and changing the number of points to check the dependency of the solution on the discretisation.

Adaptive refinement of the Gaussian process model

To identify the robust design point in a limited number of evaluations of the model m the GP-model is adaptively refined with the specific goal of identifying the most robust design point.

Therefore, the GP-model itself is used to identify regions of interest based on two criteria related to the famous compromise between exploration (low prediction confidence) and exploitation (identified areas of possible optimum). The learning function to achieve this is described in this section, starting first with an introduction of the maximum improvement function.

Maximum improvement function

The learning function introduced in this paper is based on the maximum improvement function, which was introduced in [START_REF] De Munck | An efficient response surface based optimisation method for non-deterministic harmonic and transient dynamic analysis[END_REF]. Before applying this idea to the robustness problem as stated in the previous section, the general idea is briefly summarised. The goal of the learning function is to determine which sample is the best candidate to enrich the set of calibration samples for the GP. This effectively means improving the precision of the GP around the selected sample. Here this is illustrated on a general continuous function f (u) : R nu → R, which is approximated by a GP-model g(u). Using the learning function as defined in [START_REF] Mitseas | Robust design optimization of structural systems under evolutionary stochastic seismic excitation[END_REF], a compromise is made between improving the calibration around the expected minimum using the GP mean (exploitation) and in areas of high prediction uncertainty based on the GP variance where a better minimum could be found (exploration). By iteratively enriching the calibration of the GP with the best sample improves the estimation of the minimum until a stopping criterion is eventually reached. The learning function is defined as [START_REF] De Munck | An efficient response surface based optimisation method for non-deterministic harmonic and transient dynamic analysis[END_REF]:

M I(u) = min(µ g (u)) -(µ g (u) -cσ g (u)) min(µ g (u)) , (20)
with µ g (u) the GP model prediction at u, min(µ g (u)) the current minimum, and cσ g (u) represents the variance around the prediction of u. Here, the variance is truncated at a certain confidence bound with c in Eq. ( 14). Hence, when the confidence bounds are based on, e.g., 3σ, more effort is dedicated to reducing the uncertainty about the approximation. Contrarily, lower confidence bounds, e.g., 2σ, reduce the confidence interval and favour improving approximately found maxima or minima. To identify the new candidate point u new the maximum M I is identified over the domain u ∈ U found by:

u new = argmax u∈U min(µ g (u)) -(µ g (u) -cσ g (u)) min(µ g (u)) . ( 21
)
Figure 3 shows the true function f (u) in red and the GP based approximation g(u) in black.

The black dot is a point that is part of the DOE used to calibrate the GP-model. Furthermore, this figure shows how the learning function in Eq. ( 20) is used to evaluate the point u new ∈ U to determine which point should be added to the DOE. When the GP-model is re-calibrated using the newly evaluated point u new , the minimum of f (u) is further approximated. If it is unlikely that a point u candidate provides a minimum of f (u) lower than the current min g(u), a negative improvement is obtained.

Maximum improvement of the robustness

After the introduction of the maximum improvement, the remainder of section 5 describes how this is used in this specific case of robustness. The main goal of the optimisation procedure is to identify the most robust design point in z ∈ Z, such that this design provides a minimum variation in the output interval for all x ∈ x I . This is enabled by adapting the maximum improvement, introduced in Eq. ( 20), to work directly on the minimum interval width. Specifically, it is adapted to:

M I z (z) = min z∈z I y gp (z) -y gp (z) -δ µ-σ (z) -δ µ+σ (z) min z∈z I y gp (z) -y gp (z) , (22)
y u µ gp (u) f (u) µ gp (u) + cσ gp (u) µ gp (u) -cσ gp (u) min(µ gp (u)) µ gp (u new ) -cσ gp (u new ) u new M I (unscaled)
Figure 3: Illustration of the learning function where the goal is to approximate the minimum of the true function f (u) in red, by the GP-model prediction g(u); the black dot is point that is part of the DOE and the next point u new is selected by the learning function, adapted from [START_REF] De Munck | An efficient response surface based optimisation method for non-deterministic harmonic and transient dynamic analysis[END_REF] with δ µ-σ (z)-δ µ+σ (z) the predicted minimum bound 2∆δ(z) with a confidence interval of cσ about this bound, and min z∈z I y gp (z) -y gp (z) the minimum bound predicted by the mean estimate.

Note that the mean estimated bounds correspond to 2∆y gp (z * gp ) in Figure 2b, and δ µ-σ (z)-δ µ+σ (z) to 2∆δ gp (z pot gp ) in the same figure. By reaching a M I z (z) ≤ 0, when the two intervals are equal, one can state that it is not expected with, e.g, 95% confidence for c = 1.96, that there is a smaller bound of ∆y within the current range of design parameters z ∈ Z.

Maximum improvement of the predicted bounds

The previously introduced improvement function Eq. ( 22) finds a promising design point, based on the estimates of the GP-model. However, to estimate promising design points the overall GPmodel must be refined as well, especially around these promising design points. Therefore, a second improvement function is introduced to increase the confidence of the predicted bounds.

Here the maximum improvement Eq. ( 20) is adapted to obtain a best estimate of the minimal interval width ∆δ min , which depends on both the upper and lower bound. Figure 5 illustrates the idea behind the improvement function used here. In general, the goal is to approximate the output set ỹs for each design z ∈ Z. The point that provides the largest improvement of the lower bound of this interval is given as:

M I min (θ θ θ) = min x∈x I [µ gp (θ θ θ) + cσ gp (θ θ θ)] -µ gp (θ θ θ), (23)
and the improvement of the upper bound is given as:

M I max (θ θ θ) = µ gp (θ θ θ) -max x∈x I [µ gp (θ θ θ) -cσ gp (θ θ θ)] . (24)
Note that unlike the improvement functions in Equations ( 20) and ( 22) the one given in [START_REF] Faes | Fully decoupled reliability-based design optimization of structural systems subject to uncertain loads[END_REF] and ( 24) are not normalized and calculated for each design in Z. Hence, there is a guaranteed possible improvement even if the global minimum and maximum are identified. The improvement function is illustrated in Figure 5 for a single point x * ∈ x I . In the illustrated case, the improvement of the minimum bound M I min is unlikely (negative value) while it seems likely to improve the upper limit M I max . In the end, only one candidate point can be chosen to be added to the design of experiments. Therefore, for each evaluated point the highest improvement value is used, which can either improve the lower bound or the upper bound:

M I x = max(M I min , M I max ). (25)
This means that for the illustration in Figure 5 Finally, the candidate point that performs best over the sum of the two improvement functions Eq. ( 22) and Eq. ( 25) is selected. Hence, the next candidate point θ θ θ candidate is obtained by:

θ θ θ candidate = argmax z∈Z x∈x I [M I z (z) + M I x (θ θ θ)] . (26)
Note that it is possible here to assign weighting factors to the two functions. However, to the authors knowledge no advantage is gained in this regard. Hence, these weights are not used in this work.

Stopping criterion for adaptive refinement

The role of a stopping criterion is to indicate when the algorithm reached a desired level of convergence. In this work, the stopping criterion is defined on the improvement of the robustness M I z , which means that based on the current GP-model it is unlikely to identify a point that is more robust than the current best estimate min z∈z I y gp (z) -y gp (z) . This point is identified with cσ confidence when the maximum improvement M I z ≤ 0. However, this is only achieved when the GP-model variance at location z * reduces to zero. Although possible in theory, this is highly unlikely to be achieved in practice. Hence a small error term is defined, which assures that when:

M I z ≤ , (27)
there is with 95% confidence no point R within the domain smaller than R(1+ ). Unless explicitly specified otherwise, the default value for = 1 • 10 -3 throughout this work.

Overview of the method

In Figure 6, a flowchart of the method is provided. The flowchart describes in detail the steps needed to perform the optimisation as proposed in this paper. The method starts at the initialisation where all parameters are selected by the user, i.e., correlation function, size of the initial design of experiments, value for . After this initialisation is made, the initial design of experiments is evaluated by the model m and the GP is calibrated. Hereafter, the adaptive refinement starts with finding new potential robust designs points based on the learning function in Section 5. For each newly identified point the model is evaluated m(θ θ θ candidate ) and the results are added to the Design of Experiments. This loop continues until the stopping criterion Eq. ( 27) is met. Finally, after finishing the optimisation, it is highly recommended to validate and verify the results of the GP. A good starting point to check the accuracy of the GP-model is to perform Leave-one-Out (LOO) cross-validation with the points already in the Design of experiments.

Case studies

In this section the RuLoK technique is tested and validated for different problems, which start with a set of analytical functions and build up to higher dimensional engineering examples. For the first analytical cases a comparison is made with classical optimisation techniques, which require direct evaluations of the the numerical model for each of the sample points. Moreover as the robustness measure in Eq. ( 6) requires a double-loop optimisation approach, where the outer-loop is focused on the next design point and the inner loop identifies the upper and lower bound of the response for a given design z ∈ Z. This optimisation directly uses the expensive to evaluate numerical model. Thus, the efficiency is measured in the amount of required function evaluations.

Analytical test functions

To study the basic properties of the proposed method a set of analytical test functions is used.

Each of the three test functions presents a different challenge in terms of optimisation, starting form a convex and smooth function and progressing to non-convex problems. The analytical test functions are defined as: the optimisation is to identify the value for x 1 at which the bounds on ∆f are minimal for each

f a (x 1 , x 2 ) = x 2 1 x 2 -x 2 2 , ( 28
) f b (x 1 , x 2 ) = x 2 x 1 -sin (x 1 )x 2 2 + x 2 1 , ( 29
) f c (x 1 , x 2 ) = cos(4πx 1 ) -sin(x 1 x 2 ) + x 2 , (30)
x 2 ∈ x I 2 .
This optimisation is defined as: max R(x 1 ) = min

x 1 ∈x I 1 max x 2 ∈x I 2 f n (x 1 , x 2 ) -min x 2 ∈x I 2 f n (x 1 , x 2 ) , (31)
with n indicating the three functions f a , f b , f c . In these particular cases, without the need for optimisation, one can determine that the minimum of the functions f a , f b and f c lies at

x 1 = 0; ∀x 2 ∈ x I 2 .
Nevertheless to demonstrate the additional value of the proposed method two wellknown optimisation algorithms are used in a comparison. These two optimisation approaches used in this work are: Unconstrained Optimisation (UO) where the minimum of a function is searched using a quasi-Newton algorithm; an other strategy is to use a Generic Algorithm (GA) to solve the outer-loop where the bounds of the response in the inner-loop are identified using UO. The population for the GA is set to a default value of 20.

The results of the method and these of the classical optimisation approaches are compared in Table 1. It is noticed that the proposed method outperforms the brute optimisation approaches, which is expected with the use of a meta-model. The table also shows that depending on the level of confidence the number of iterations increases. Note that the amount of iterations needed to obtain a result is difficult to estimate a priori as this depends on the underlying problem and the correctness of the GP model at each iteration. The error term in the table refers to the discretization error introduced by using a fixed grid to sample the meta model. For both function f a and f b the optimal point is part of the samples in the grid using n samples = 501. However, for function f c this dependence is checked and the optimal point is not part of the grid points n samples = 200 or n samples = 500. Therefore, the analysis returns the next best point, which is the closest to the optimal point. Using a larger number of grid-points will therefore increase the accuracy of the estimation at a higher computational cost. To further illustrate how the method works Figure 7 shows the function value for all three functions f a ,f b and f c at each design point x 1 . For each function the true bounds are given by the black dashed lines, the evaluated point are indicated with a green cross, and the predicted upper-and lower-bound are given in red and blue, including their 95% confidence intervals, and the optimal design point is indicated by a circle. Starting at the top of Figure 7 Finally, the bottom graph of Figure 7 shows the more complex function f c with the optimum at x 1 = 0. This figure illustrates the additional function evaluations needed to ensure the global minimum was found, and not one of the many local minimums. Note that for this case the optimal robust point x 1 = 0 is not part of the grid as the grid is discretized by an even number of samples, which include the end and start point.

Function Method Optimum Iterations Evaluations Error

* confidence f a analytic 0 - - - - f a RuLoK 0 36 38 ≈ 0.02 1, 96σ f a RuLoK 0 42 44 ≈ 0.02 3σ f a UO -7e-6 3 184 - - f b analytic 0 - - - - f b RuLoK 0 28 30 ≈ 0.02 1, 96σ f b RuLoK 0 49 51 ≈ 0.02 3σ f b GA 0.099 77 50680 - - f c analytic 0 - - - - f c RuLoK -2.
In Figure 8 the meta-model of function f a is plotted with the black surface the mean response, the red and blue surfaces the lower-and upper-bound of the 95% confidence intervals, and the green dots are the points used to calibrate the GP-model. This figure illustrates the dispersion of the evaluation points at the edges of the domain and concentration of points around the optimal point, which reduces the variance of the GP-model is this location. Hence, the distance between the bounds increases in locations that are further from the optimal point since there are considerably less points evaluated here. Nevertheless, it is possible to use the GP-model further to analyse the problem at hand. However, one should be aware that due to the selection of training points an overall agreement between the GP-model and the underlying problem is not guaranteed.

Plate subjected to a point load

In this case study, the thickness of a plate with two equal sides of 100mm is chosen within the interval t ∈ [START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF][START_REF] Moore | Interval analysis[END_REF]mm. The uncertain parameter is the Young's modulus of the material, which At the start of the analysis two initial points are evaluated based on Latin Hyper-Cube sampling. Figure 9 shows the results that are obtained after just 9 function evaluations, with the true bounds of the model in dashed black lines, the GP-model prediction of the upper-and lowerbound in red and blue, and the confidence intervals as colored areas. The optimal design point t robust = 6mm is as expected, the thickest plate. The rationale behind this simple example is that the thickest plate will bend less than a thinner plate under identical uncertainty of the Young's modulus. However, Figure 9 presents an illustration of the refinement around the optimal point, with only two evaluation points lower than 4.5mm. The order of the points that are added is further highlighted by the numbers next to the crosses in the plot starting with the initial evaluations 1 and 2, up to 9, the final point.

Although the physical interpretation of the problem explains the identified optimum a doubleloop approach is used to validate this result. Here using CO a total of 82 evaluations of the numerical model were needed to identify the optimum t robust = 6mm, which is identical. However, with this classical optimisation no additional information is obtained regarding the problem that is studied.

The borehole function

The second engineering example is the Borehole function [START_REF] Surjanovic | Virtual library of simulation experiments: Test functions and datasets[END_REF], which is a typical test case for computer experiments. The borehole function describes the water flow f borehole though a borehole between two underground aquifers by the flow rate of the water m 3 /year: 2. Two cases are considered with this example, first a case where only two parameters are uncertain and the others are taken at the midpoint, second a case where all parameters are considered uncertain. The physical interpretation of the location of the robust point is that a borehole with a smaller diameter limits the possible flow through the borehole. However, for the length of the borehole this observation is not obvious. The results in Figure 10 are obtained with a total of 35 evaluations of the borehole function including the four initial evaluations.

Borehole function with six uncertain parameters

In this case all six uncertain-and the two design-parameters of the previously discussed borehole function Eq. ( 32) are considered within the ranges as defined in Table 2. The results of the analysis are shown in Figure 11, which shows the true interval width on the top and the GP prediction on the bottom left and the interval width based on the 95% CI on the right. The number of evaluations to obtain these results has only increased slightly to 64, which includes 8 initial evaluations, while the complexity of the problem is increased by four additional uncertain parameters. The location of the robust design point remained at the lower-bound of the diameter and the upper-bound of the borehole length. The physical reason for this difference is not directly clear from the formulation of the borehole function. However, the additional parameters seem to have little effect to the overall behaviour of the function while the width of the interval has increased slightly, which can be seen by comparing Figure 10 and Figure 11. To better understand the effect of the additional parameters the interval sensitivities are investigated. The reader is referred to [START_REF] Moens | Interval sensitivity theory and its application to frequency response envelope analysis of uncertain structures[END_REF] for a thorough discussion about interval sensitivities. However, note that the fundamental difference between the classical sensitivity studies and interval sensitivities is that the latter is valid over the full range of the interval, while the former focuses on local sensitivities, which are not valid over the full range of the interval. The interval sensitivities for the borehole function with six uncertain parameters are provided in Figure 12, which shows that the radius of influence r, transmissivity of the upper aquifer T u , and lower aquifer T l have an negligible effect on the output interval. Moreover, this figure shown that all parameters behave the least sensitive around the robust design point. The latter means that with a relative change of input interval width only a minimal change in output interval width happens. Although the obtained results are convincing and could be compared with the true solution, this is not always possible especially with the use of complex numerical models. However, one can validate the GP-model based on the points that were evaluated in the Design-of-Experiments, which provides an indication about the correctness to capture the underlying physical behaviour.

This validation is accomplished by a number of tests shown in Figure 13 which are based on the Leave-One-Out prediction of the points within the DOE. Note that this is a conservative choice We can conclude that the GP-model performs well at low flow rates, with an error that increases at higher flow rates.

Discussion

In general the results of the presented method are convincing and show that this method is capable of identifying the robust design point with only a limit amount of evaluations of the underlying expensive function, which is demonstrated in a number of case studies. Nevertheless, a few things are noted by the authors that should be addressed for further research and implementation. As mentioned before, the obtained results are based on the GP model as implemented in UQlab [START_REF] Marelli | UQLab: A Framework for Uncertainty Quantification in Matlab[END_REF] for all case studies. It is noted by the authors that using different implementations of the GP can lead to an increase in the number of iterations before convergence is reached. This is attributed to the use of a noise parameter in the GP, which is set at a minimum of 1e -4 for the Matlab build-in implementation [START_REF]Matlab and statistics and machine learning toolbox release[END_REF]. Hence, the error term in Equation 27 should increase to reflect this.

The number of samples in the initial DOE can effect the convergence and in this paper, as a rule of thumb, the number of initial evaluations is kept at the total amount of uncertain-and design-parameters. Quantifying the effect of the initial population size on the rate on convergence is challenging as this depends on the underlying problem, i.e., that what is resembled by the GP model. This rule of thumb is regarded as the minimal amount of initial evaluations needed by the GP to make a first estimation. Nevertheless, the number of iterations is difficult to determine a priory, as this depends on the complexity of the response surface, the added value of the point added at each iteration, and the calibration error of the GP model.

Finally, as the improvement function is evaluated on a fixed number of equally spaced grid points, a limited precision is reached. Although using a large number of grid points the precision increased, the computational cost to evaluate all these points increases exponential in d-dimensions O(n -d ) for a full grid. Hence, in high dimensional cases this becomes a bottleneck without sacrificing the resolution of the grid.

Conclusion

This paper introduces a novel method to design robust structures in an early stage of development under lack-of-knowledge uncertainty. The presented method uses an adaptively refined GP-model to perform the global optimisation of the robustness and locate the most promising designs, which are the least sensitive to the modelled sources of uncertainty. Based on a set of analytical test functions the effectiveness and efficiency of the proposed method is demonstrated and compared with typical well-known optimisation algorithms. It is shown that the proposed Future research is aimed at enlarging the application domain of the proposed method, specifically for time-dependant problems.

  a) Illustration of the Gaussian Prosses model spanning the uncertain x, design z and output y space, with the GP mean prediction µgp in black and the µ ± cσ CI in red and blue surfaces; Two potential design points z * gp and z pot gp are shown as a slice.

  b) Illustration of the (z)-(y) perspective where the Confidence bounds around the upper-bound y gp (z) and lower-bound y gp (z) are given by the red and blue area's; and the two designs z * gp and z pot gp show the difference in mean bound prediction ∆ygp(z * gp )and the lowest CI prediction ∆δgp(z pot gp ) of the bounds

Figure 2 :

 2 Figure 2: Illustration of the domain to determine the robustness based on the GP-model predictions

Figure 4 Figure 4 :

 44 Figure 4: Illustration of the predicted mean bound ∆y gp (z) = y gp (z) -y gp (z) and the minimum bound based on the confidence interval ∆δ(z) = δ µ+σ (z) -δ µ-σ (z)

  only the value of M I max is saved for the point x * . y x M I max (z, x * ) M I min (z, x * ) x * ∆y z ∆δ min ∆δ max y gp y gp + cσ gp y gp -cσ gp

Figure 5 :

 5 Figure 5: Illustration of the learning function for a candidate point x * , showing the MI of the lower and upper bound; here the improvement of the lower bound is negative

Figure 6 :

 6 Figure 6: Flowchart of the robustness under lack-of-knowledge method

  function f a is shown where the gradient decreases when moving towards the robust design point x 1 = 0. The middle sub-figure illustrates the function f b with larger confidence bounds around the predicted optimum, shown by the red and blue areas. It is also shown that the confidence interval about the upper bound is larger than this of the lower bound, which is exactly the goal during optimisation.

Figure 7 :

 7 Figure 7: The GP predicted bounds of the interval valued uncertainty including the 95% confidence intervals; for from the top to the bottom function f a ,f b and f c , respectively.

Figure 8 :

 8 Figure 8: GP-model prediction of function A with the black surface the mean response, the red and blue surfaces the 95% CI on the mean prediction, and the green points indicating the evaluated points

Figure 9 :

 9 Figure 9: GP predicted bounds of the interval valued uncertainty including the order in which the points within the DOE where evaluated; including the 95% confidence intervals, evaluated points indicated by crosses, and the optimal design point indicated by a circle.

6. 3 . 1 .

 31 Borehole function with two uncertain parametersIn this first case only the potentiometric head of the upper aquifer H u and the hydraulic conductivity K w are regarded as uncertain. The remaining uncertain parameters are taken at the midpoint of their interval. The results of the analysis are shown in Figure10, which shows a contour plot of the true interval width on the top, the predicted interval width based on the mean of the GP-model below, and the minimal interval width based on the 95% CI next to it.In all contour plots of Figure10the red circle and green dot indicate the location of the robust design point, located at the lower-bound of the diameter and the upper-limit of the length of the borehole. In addition the blue dots indicate the points where the original function was sampled.

Figure 12 :

 12 Figure 12: Relative interval sensitivity of the uncertain parameters r, T u , H u , T l , H l and K w in function of the borehole diameter r w and length L.

Figure 13 :

 13 Figure 13: Validation tests for the borehole function with two uncertain parameters: (a) shows the cross-validated prediction vs. the true function value, (b) shows the standard normalized residuals of the cross-validated GP model within the 95% bounds in red, and (c) shows the cross validated quantiles vs. the true quantiles

Table 1 :

 1 Results of the analytic test functions

			5e-2	279	281	0.05	1, 96σ
	f c	RuLoK	1e-2	172	174	0.02	1, 96σ
	f c	RuLoK	0	214	216	≈ 0.02	1, 96σ
	f c	RuLoK	-1e-2	242	244	0.02	3σ
	f c	GA	4.4e-5	30	2760857	-	-

* the discretization error of the grid is determined by ∆x/n points = 10/501 for f a and f b .

  Figure 10: top: contour plot of the true interval width in function of the design parameters, borehole radius r w and the borehole length L and only the potentiometric head of the upper aquifer H u and hydraulic conductivity K w

	2 1	3 4	4 6	5 9	7 1 8 3 9 6
	2 1	3 4	4 6	5 9	7 1	8 3	9 6	1 0 8		2 1	3 4	4 6	5 9	7 1 8 3 9 6
	are regarded uncertain; bottom: mean GP prediction of the interval width (left) and the minimal interval width
	by 95% CI (right)
	parameter										x		x	x	unit
	radius of influence									r I 100		50 000 2550	m
	transmissivity of upper aquifer						T I u	63 070 50 000 56 535 m 2 /year
	potentiometric head of the upper aquifer	H I u	990		1110	1050	m
	transmissivity of lower aquifer						T I l	63.1		116	89.55	m 2 /year
	the potentiometric head of the lower aquifer H I l	700		820	760	m
	hydraulic conductivity of the borehole	K I w	9855		12 045 10 950 m 2 /year

Table 2 :

 2 Parameters of the borehole function

  Figure11: top: contour plot of the true interval width in function of the design parameters: borehole radius r w , and the borehole length L; with uncertain parameters: the radius of influence r, transmissivity of upper aquifer T u , the potentiometric head of the upper aquifer H u , transmissivity of lower aquifer T l , the potentiometric head of the lower aquifer H l , the hydraulic conductivity of the borehole K w ; bottom: mean GP prediction of the interval width (left) and the minimal interval width by 95% CI (right)

	3 6	5 8	7 9	1 0 0	1 2 1	1 4 2 1 6 4
	3 6	5 8	7 9	1 0 0	1 2 1	1 4 2	1 6 4	3 6	5 8	7 9	1 0 0	1 2 1	1 4 2	1 6 4

Acknowledgements

The authors gratefully acknowledge the support of the Research Foundation Flanders (FWO) under grant 1SA3919N (C. van Mierlo) and GOC2218N (A. Persoons).


Fabien Soulier

Fanny Le Floch

Serge Bernard

Guy Cathébras

New Dependability Approach for Implanted Medical Devices

Functional Electrical Stimulation (FES) is an attractive solution to restore some lost or failing physiological functions. Obviously, the FES system may be hazardous for patient and the reliability and dependability of the system must be maximal. Unfortunately, the present context, where the associated systems are more and more complex and their development needs very cross-disciplinary experts, is not favorable to safety. Moreover, the direct adaptation of the existing dependability techniques from domains such as space or automotive is not suitable.

Firstly, this paper proposes a strategy for risk management at system level for FES medical implant. The idea is to give a uniform framework where all possible hazards are highlighted and associated consequences are minimized. Then, the paper focuses on one of the most critical part of the FES system: analog micro-circuit which generates the electrical signal to electrode. As this micro-circuit is the closest to the human tissue, any failure might involve very critical consequences for the patient. We propose a concurrent top-down and bottom-up approach where the critical elements are highlighted and an extended risk analysis is performed.

I. INTRODUCTION

Functional Electrical Stimulation (FES) is quite often used to restore damaged physiological functions. Because of the implantation of the system inside the human body, the reliability and safety of the implant have to be perfectly controlled. Unfortunately, the complexity of implant endlessly increases, as a result, providing the highest safety is a real challenge. As an illustration, until now, the classical cochlear implant was composed of an external part including a microphone, a speech processor and an RF transducer coupled with a coil implanted beneath the skull's skin and an internal circuit inside the cochlea stimulating auditory nerves with an electrical current. The system was rather basic and the dependability study was easily done. But with the increase of the complexity and versatility of the system, the dependability study [START_REF] Halperin | Security and privacy for implantable medical devices[END_REF]- [START_REF] Knight | Safety critical systems: challenges and directions[END_REF] has to be taken more seriously.

The development of such an implant brings a very highly cross-disciplinary context where physicians, surgeons, microelectronics and mechanical specialists have to work together. They have to take into account the expectation of each other even if the communication between these different communities is not easy. For example, the surgeon needs to be educated on how to manipulate an implant, which is highly sensitive to electric discharge. The experience shows that sometimes working implant before the surgery was defective once they have been embedded inside the patient, due to a bad manipulation. The solution is then to test the implant before ending the surgery.

Firstly, this paper presents a new strategy for risk management at system level for FES medical implant. The aim is to propose a uniform framework to highlight each potential hazard and solution to control and prevent them. The existing techniques developed in aeronautic, space or even automotive fields are not totally appropriated for the medical domain. The third section focuses on the electronic part of the cochlear implant. Finally, the last section introduces an efficient concurrent top-down and bottom-up approach on the analog part of the circuit.

II. STRATEGY FOR RISK MANAGEMENT AT SYSTEM LEVEL

A. Objectives

The high-level dependability management [START_REF] Villemeur | sûreté de fonctionnement des systèmes industriels[END_REF] gives us a chance to consider every potential defective point. One of the main objectives is to give a set of tools and guide lines suitable for any expert or set of experts of the entire system. The base of the proposed risk management is an algorithm which described the successive steps. This algorithm must be as simple as possible to achieve the best coverage of defective cases and propose the highest level of confident and safety of the system.

B. Algorithm for risk management

As shown in figure 1, we have developed an iterative algorithm for the dependability management. First, the functional analysis defines the system, its environment and its external limits. It points up that not only the implant inside the human body has to be taken care, but also the patient itself, the surgeon, the practitioner, the patient's activities. . . It is useful to understand and determine the patient's way of life. What sport does he practice? (Shocks can damage the implant) Does he need specific settings if he works in a noisy surrounding?

This step is essential because if the definition of the expected system functionalities is not precise and the possible user cases and external constraints are not exhaustive, a potential hazard context might appeared because it was not taken into account. The other objective of this first step is to break up the system into functional blocks and general properties. Theses subsystems allow us to define working groups associated with specific expertise.

The second step is the hazard's analysis which draws up the full list of hazards that can occur. By means of a chart, we collect hazards into two separates categories: a specific and a generic part. The specific part corresponds with the functional decomposition of the project: the stimulation architecture, the wireless and power interface, the embedded software, the neuro-sensing interface, the micro-mechanics module and the battery. The generic part is the segmentation of the system into properties such as: chemical properties, mechanical, physiological (ex: biocompatibility of the cover material of the implant with the human tissue), with regards to the environment (ex: electromagnetic compatibility), the diagnosis, concerning the radiation protection, the informations given by the manufacturer, concerning the practitioner, the surgeon, the patient, the implantation, the implant-life inside the human body [START_REF] Creasy | Implications of implantation of bladder stimulation systems in children and adolescents[END_REF], the disposal procedures. Once the patient is dead, the implant is removed from the body for environmental reasons. This preliminary version is written by the person in charge of the dependability study. Once it is established, it is then sent to the working team, separated in smaller teams depending on their skills. The hazard chart has to be as exhaustive as possible. Afterwards, a risk analysis will be done: the purpose is to estimate the dependability level before and after the countermeasure to check the efficiency of the taken modification. For instance, there is a level shifter (see figure 3) which controls the high voltage part of the circuit, thanks to a low voltage digital part. We found out that this was a sensitive part of the implant, so the countermeasure we are developing is a simple down shifter. If the comparison of its input and the output of the level shifter is different, a warning is sent to the logic control, and the implant switches in a fail-soft mode.

The dependability study has to carry on over the entire life of the project, thanks to the experience feedback, the study can be completed and improved. The link between the person in charge of the dependability study and the after sales services is essential to keep collecting information for the existing implant and the future generation of implant.

III. RISK ANALYSIS OF THE ELECTRONIC PART OF THE

IMPLANT

In this section, we focus on one of the most critical part of the implant: the electronic part. The analog part of the microcircuit generates the electrical signal send to the electrodes. As this micro-circuit is the closest to the human tissue, any failure might involve very critical consequences for the patient. As an illustration, we apply our approach to the cochlear implant.

A. Case study: the cochlear implant

Cochlear implant is used for deafness. The fundamental principle is to convert sound recorded with a microphone to electrical stimulations inside the cochlea. In this section, we present an overview of the future generation of this kind of implant where, to make life easier for the patient, the entire system will be implanted into the human body. The device is composed of an external controller [START_REF] Cherukuri | Biosec: A biometric based approach for securing communication in wireless networks of biosensors implanted in the human body[END_REF], which has to be used every two days to recharge the internal integrated circuit by means of wireless link, as shown on figure 2. The major difficulty with the battery is to give to the patient the freedom to recharge the implant whenever he wants (as long as he does it at least once a week), instead of restraining him on a precise time.

The embedded part inside the body is divided into two subsystems: a digital and an analog module. The digital module manages the communication with the external controller. It also receives the sound picked up from a microphone implanted inside the cochlear. The sound is processed before being send to the logic control, which determines the shape, the amplitude and the duration of the stimulation current [START_REF] Guiraud | An implantable neuroprosthesis for standing and walking in paraplegia: 5-year patient follow-up[END_REF], [START_REF] Guiraud | Original electronic design to perform epimysial and neural stimulation in paraplegia[END_REF].

The power is supplied by a battery, implanted beneath the collarbone. The analog module contains a digital to analog converter (DAC), which provides the stimulation current. The high voltage generator is necessary to reach the maximum

B. The analog part of the FES implant

Having a failure from the analog part of the micro-circuit can be a disaster because among other things, it is in direct contact with the human tissue. We use a approach based on concurrent top-down and bottom-up methods (see figure 4). The tangle of these two approaches brings out every potential failure and give us the exhaustive list of failures and allows us to develop specific Built-In-Self-Test and Built-In-Self-Repair solutions to enhance the safety during the life of the system.

1) The top-down approach: The top-down approach consists of breaking down the system into finer details. We refine the system until we get to elementary structures.(The division stops when it seems irrelevant to go further). Each step of decomposition is analyzed in terms of failures, so the advantage is to consider every potential failure.

The study starts once the functional description has been done. The global view gives the opportunity to work not only on the embedded system inside the patient, but also, on every human which interact with the implant (the patient, the surgeon, the practitioner,. . . ). The implant life is also examined: it begins with its manufacture, carriage, storage, implantation, life in the human body, and finally, its disposal procedures. Then, each specific expert team looks downer inside the global system. This approach helps us to cover the system from a macro overview to a micro decomposition.

2) The bottom-up approach: The bottom-up approach starts with a basic element of the system, known as deficient and the following consequences of the deficiency are tracked up to the actual hazard. We can apply this method because the analog circuit is an improvement of an older version, so we can use the experience feedback. As an example, the impedance measurement helps to check if the electrodes are in open circuit, or short circuit, or even damaged. This is a countermeasure that has been implemented in the latest version of the implant, we realized that if a short circuit occurred, we were not able to detect that the failure came from the electrode, we had to change not only the electrode, but the entire implant. Replacing a defective implant is not a minor operation. The patient has to go through another surgery and there is a delay to respect between the explantation of the defective implant and the implantation of the new one. Once every critical point has been listed, we can treat them by applying efficient countermeasure.

IV. CONCLUSION

The complexity and versatility of implant leads to develop new methods proper to the medical field. One of the main tools for safety improvement is the algorithm of the risk management which gives a guideline to follow in order to obtain an exhaustive list of potential failures. Once this list is established, we can develop efficient Built-In-Self-Test and Built-In-Self-Repair solutions. The control of every critical element of the implant and its environment helps to increase the reliability and safety of the system. This approach is iterative, and needs to be updated during the entire life of the implant. In the next future, we will focus our research on the analog part of the circuit, which has been technically improved. Thanks to the experiment feedback on the old implant generations and the developed algorithm, this coming and going approach will provide a complete list of failures. Indeed, the top-down approach covers the system from a macro overview to a micro decomposition, whereas the bottom-up approach uses the experience feedback to repair know defects. This merge will bring us a highly safe and reliable medical implant to the patient advantage.

Fig. 1 .

 1 Fig. 1. Algorithm of the dependability system management.

Fig. 2 .

 2 Fig. 2. Overview of the implant.

Fig. 3 .

 3 Fig. 3. Schematic of level shifter.

Fig. 4 .

 4 Fig. 4. Top-down and bottom-up method.