In engineering design, it is important to guarantee that the values of certain quantities such as stress level, noise level, and vibration level, stay below a certain threshold in all possible situations, i.e., for all possible combinations of the corresponding internal and external parameters. Usually, the number of possible combinations is so large that it is not possible to physically test the system for all these combinations. Instead, we form a computer model of the system and test this model. In this testing, we need to take into account that the computer models are usually approximate. In this paper, we show that the existing techniques for taking model uncertainty into account overestimate the uncertainty of the results. We also show how we can get more accurate estimates.

## Introduction

### Bounds on Unwanted Processes: An Important Part of Engineering Specifications.

An engineering system is designed to perform certain tasks. In the process of performing these tasks, the system also generates some undesirable side effects: it can generate noise, vibration, heat, and stress.

We cannot completely eliminate these undesired effects, but specifications for an engineering system usually require that the size $q$ of each of these effects does not exceed a certain predefined threshold (bound) $t$. It is therefore important to check that this specification is always satisfied, i.e., $q\u2264t$ in all possible situations.

### How Can We Check That Specifications are Satisfied for All Possible Situations: Simulations are Needed.

To fully describe each situation, we need to know which parameters $p1,p2,\u2026$ characterize this situation, and we need to know the exact values of all these parameters. These may be external parameters such as wind speed and load for a bridge. This may be internal parameters such as the exact value of the Young’s module for a material used in the design.

Thus, without losing generality, we can always assume that the set of possible values of each parameter $pi$ is given by Eq. (1).

We would like to make sure that the quantity $q$ satisfies the desired inequality $q\u2264t$ for *all* possible combinations of values $pi\u2208[p_i,p\xafi]$. Usually, there are many such parameters, and thus, there are many possible combinations—even if we limit ourselves to extreme cases, when each parameter $pi$ is equal to either $p_i$ or $p\xafi$, we will still get $2n$ possible combinations. It is therefore not feasible to physically check how the system behaves under all such combination. Instead, we need to rely on computer simulations.

### Formulation of the Problem.

There are known techniques for using computer simulation to check that the system satisfies the given specifications for all possible combinations of these parameters; see, e.g., [5] and references therein. These techniques, however, have been originally designed for the case in which we have an exact model of the system.

In principle, we can also use these techniques in more realistic situations, when the corresponding model is only approximate. However, as we show in this paper, the use of these techniques leads to overestimation of the corresponding uncertainty. We also show that a proper modification of these techniques leads to a drastic decrease of this overestimation and thus to more accurate estimations.

## How to Check Specifications When We Have an Exact Model of a System: Reminder

### Case of an Exact Model: Description.

### In Most Engineering Situations, Deviations From Nominal Values are Small.

Usually, possible deviations $\Delta pi=defpi\u2212p\u02dci$ from nominal values are reasonably small; see, e.g., [6]. In this paper, we will restrict ourselves to such situations.

### How to Use the Linearized Model to Check That Specifications are Satisfied: Analysis of the Problem.

To make sure that we always have $q\u2264t$, we need to guarantee that the largest possible value $q\xaf$ of the function $q$ does not exceed $t$.

### How to Estimate the Derivatives ci

In many practical cases, we have an explicit equation or, more generally, a known program for computing $q(p1,\u2026,pn)$. In this case, by explicitly differentiating the corresponding expression—or by applying an automatic differentiation algorithm to the corresponding program—we can get equations for computing the derivatives $ci$.

In many real-life situations, however, in our computations, we use proprietary software—for which the corresponding program is not available to us. In such situations, we cannot use automatic differentiation tools, and we can use only the results of applying the algorithm $q$ to different tuples.

Thus, we arrive at the following technique (see, e.g., [7]).

### How to Use the Linearized Model to Check That Specifications are Satisfied: Resulting Technique.

We know:

- •
An algorithm $q(p1,\u2026,pn)$ that, given the values of the parameters $p1,\u2026,pn$, computes the value of the quantity $q$.

- •
A threshold $t$ that needs to be satisfied.

- •
For each parameter $pi$, we know its nominal value $p\u02dci$ and the largest possible deviation $\Delta i$ from this nominal value.

Based on this information, we need to check whether $q(p1,\u2026,pn)\u2264t$ for all possible combinations of values $pi$ from the corresponding intervals $[p\u02dci\u2212\Delta i,p\u02dci+\Delta i]$.

We can perform this checking as follows:

- 1.
First, we apply the algorithm $q$ to compute the value $q\u02dc=q(p\u02dc1,\u2026,p\u02dcn)$;

- 2.Then, for each $i$ from 1 to $n$, we apply the algorithm $q$ to compute the value$qi=q(p\u02dc1,\u2026,p\u02dci\u22121,p\u02dci+\Delta i,p\u02dci+1,\u2026,p\u02dcn);$
- 3.
After that, we compute $q\xaf=q\u02dc+\u2211i=1n|qi\u2212q\u02dc|$; and

- 4.
Finally, we check whether $q\xaf\u2264t$.

If $q\xaf\u2264t$, this means that the desired specifications are always satisfied. If $q\xaf>t$, this means that for some combinations of possible values $pi$, the specifications are not satisfied.

### Possibility of a Further Speedup.

Equation (9) requires $n+1$ calls to the program that computes $q$ for given values of parameters. In many practical situations, the program $q$ takes a reasonably long time to compute, and the number of parameters is large. In such situations, the corresponding computations require a very long time.

The possibility to use Cauchy distributions comes from the fact that they have the following property: if $\eta i$ are independent variables that are Cauchy distributed with parameters $\Delta i$, then for each tuple of real numbers $c1,\u2026,cn$, the linear combination $\u2211i=1nci\xb7\eta i$ is also Cauchy distributed, with parameter $\Delta =\u2211i=1n|ci|\xb7\Delta i$.

Thus, we can find $\Delta $ as follows [8]:

- 1.
First, for $k=1,\u2026,N$, we simulate random variables $\eta i(k)$ which are Cauchy distributed with parameters $\Delta i$.

- 2.For each $k$, we then estimate $\Delta q(k)=\u2211i=1nci\xb7\eta i(k)$ as $\Delta q(k)=q(k)\u2212q\u02dc$, where(12)$q(k)=q(p\u02dc1+\eta 1(k),\u2026,p\u02dcn+\eta n(k))$
- 3.Based on the population of $N$ values $\Delta q(1),\u2026,\Delta q(N)$ which are Cauchy distributed with parameter $\Delta $, we find this parameter, e.g., we can use the maximum-likelihood method according to which the desired value $\Delta $ can be found from the following equation:which can be easily solved by bisection if we start with the interval $[\Delta _,\Delta \xaf]$ in which $\Delta _=0$ and $\Delta \xaf=maxk|\Delta q(k)|$;$\u2211k=1N11+(\Delta q(k))2\Delta 2=N2$
- 4.
Finally, we follow Eq. (6) and compute $q\xaf=q\u02dc+\Delta $ (see, e.g., [8] for technical details).

In this Monte-Carlo-type technique, we need $N+1$ calls to the program that computes $q$. The accuracy of the resulting estimate depends only on the sample size $N$ and *not* on the number of inputs $n$. Thus, for a fixed desired accuracy, when $n$ is sufficiently large, this method requires much fewer calls to $q$ and is, thus, much faster. For example, if we want to estimate $\Delta $ with relative accuracy 20%, then we need $N=100$ simulations, so for $n\u226b200$, this method is much faster than a straightforward application of Eq. (9).

### For Many Practical Problems, We Can Achieve an Even Faster Speedup.

In both methods described in this section, we apply the original algorithm $q(p1,\u2026,pn)$ several times: first, to the tuple of nominal values $(p\u02dc1,\u2026,p\u02dcn)$ and then, to several other tuples $(p\u02dc1+\eta 1,\u2026,p\u02dcn+\eta n)$. For example, in the linearized method (9), we apply the algorithm $q$ to tuples $(p\u02dc1,\u2026,p\u02dci\u22121,p\u02dci+\Delta i,p\u02dci+1,\u2026,p\u02dcn)$ corresponding to $i=1,\u2026,n$.

*linear*equations

This idea—known as *local sensitivity analysis*—is successfully used in many practical applications; see, e.g., [9,10].

*Comment.* As we have mentioned earlier, in this paper, we consider only situations in which the deviations $\Delta pi$ from the nominal values are small. In some practical situations, some of these deviations are not small. In such situations, we can no longer use linearization, and we need to use global optimization techniques of *global sensitivity*; see, e.g., [9,10].

## What If We Take Into Account Model Inaccuracy

### Models are Rarely Exact.

Engineering systems are usually complex. As a result, it is rarely possible to find explicit expressions for $q$ as a function of the parameters $p1,\u2026,pn$. Usually, we have some approximate computations. For example, if $q$ is obtained by solving a system of partial differential equations, we use, e.g., the finite element method to find the approximate solution and thus, the approximate value of the quantity $q$.

### How Model Inaccuracy is Usually Described.

### How This Model Inaccuracy Affects the Above Checking Algorithms: Analysis of the Problem.

Hence, we arrive at the following method.

### How This Model Inaccuracy Affects the Above Checking Algorithms: Resulting Method.

We know:

- •
An algorithm $Q(p1,\u2026,pn)$ that, given the values of the parameters $p1,\u2026,pn$, computes the value of the quantity $q$ with a known accuracy $\u03f5$.

- •
A threshold $t$ that needs to be satisfied.

- •
For each parameter $pi$, we know its nominal value $p\u02dci$ and the largest possible deviation $\Delta i$ from this nominal value.

Based on this information, we need to check whether $q(p1,\u2026,pn)\u2264t$ for all possible combinations of values $pi$ from the corresponding intervals $[p\u02dci\u2212\Delta i,p\u02dci+\Delta i]$.

We can perform this checking as follows:

- 1.First, we apply the algorithm $Q$ to compute the value(27)$Q\u02dc=Q(p\u02dc1,\u2026,p\u02dcn)$
- 2.Then, for each $i$ from 1 to $n$, we apply the algorithm $Q$ to compute the value(28)$Qi=Q(p\u02dc1,\u2026,p\u02dci\u22121,p\u02dci+\Delta i,p\u02dci+1,\u2026,p\u02dcn)$
- 3.After that, we compute(29)$B=Q\u02dc+\u2211i=1n|Qi\u2212Q\u02dc|+(2n+1)\xb7\u03f5$
- 4.
Finally, we check whether $B\u2264t$.

If $B\u2264t$, this means that the desired specifications are always satisfied. If $B>t$, this means that we cannot guarantee that the specifications are always satisfied.

*Comment 1.* Please note that, in contrast to the case of the exact model, if $B>t$, this does not necessarily mean that the specifications are not satisfied: maybe they are satisfied, but we cannot check that because we know only approximate values of $q$.

*Comment 2.* Similar bounds can be found for the estimates based on the Cauchy distribution.

*Comment 3.* The above estimate $B$ is not the best that we can get, but it has been proven that computing the best estimate would require unrealistic exponential time [11,12]—i.e., time that grows as $2s$ with the size $s$ of the input; thus, when we consider only feasible algorithms, overestimation is inevitable.

*Comment 4.* Similar to the methods described in the previous section, instead of directly applying the algorithm $Q$ to the modified tuples, we can, wherever appropriate, use the above-mentioned local sensitivity analysis technique.

*Problem.* When $n$ is large, then, even for reasonably small inaccuracy $\u03f5$, the value $(2n+1)\xb7\u03f5$ is large.

In this paper, we show how we can get better estimates for the difference between the desired bound $q\u02dc$ and the computed bound $Q\xaf$.

## How to Get Better Estimates

### Main Idea.

Model inaccuracy comes from the fact that we are using approximate methods to solve the corresponding equations.

Strictly speaking, the resulting inaccuracy is deterministic. However, in most cases, for all practical purposes, this inaccuracy can be viewed as random: when we select a different combination of parameters, we get an unrelated value of inaccuracy.

### Technical Details.

What is a probability distribution for these random variables?

All we know about each of these variables is that its values are located somewhere in the interval $[\u2212\u03f5,\u03f5]$. We do not have any reason to assume that some values from this interval are more probable than others, so it is reasonable to assume that all the values are equally probable, i.e., that we have a uniform distribution on this interval.

For this uniform distribution, the mean is 0, and the standard deviation is $\sigma =\u03f5/3$.

### Auxiliary Idea: How to Get a Better Estimate for q˜.

The left-hand side is the arithmetic average of $n+2$ independent identically distributed random variables, with mean 0 and variance $\sigma 2=\u03f52/3$. Hence (see, e.g., [13]), the mean of this average $\Delta q\u02dcnew$ is the average of the means, i.e., 0, and the variance is equal to $\sigma 2=\u03f52/[3\xb7(n+2)]\u226a\u03f52/3=\sigma 2[\Delta q\u02dc]$.

Thus, this average $Q\u02dcnew$ is a more accurate estimation of the quantity $q\u02dc$ than $Q\u02dc$.

### Let us Use This Better Estimate for q˜ When Estimating the Upper Bound q¯.

Let us estimate the accuracy of this new approximation.

Here, the differences $\Delta q\u02dcnew$ and $\Delta qi$ are independent random variables. According to the central limit theorem (see, e.g., [13]), for large $n$, the distribution of a linear combination of many independent random variables is close to Gaussian. The mean of the resulting distribution is the linear combination of the means, thus equal to zero.

Here, inaccuracy grows as $2n+1$, which is much better than in the traditional approach, where it grows proportionally to $2n+1$—and we achieve this drastic reduction of the overestimation, basically by using one more run of the program $Q$ in addition to the previously used $n+1$ runs.

So, we arrive at the following method.

### Resulting Method.

We know:

- •
An algorithm $Q(p1,\u2026,pn)$ that, given the values of the parameters $p1,\u2026,pn$, computes the value of the quantity $q$ with a known accuracy $\u03f5$.

- •
A threshold $t$ that needs to be satisfied.

- •
For each parameter $pi$, we know its nominal value $p\u02dci$ and the largest possible deviation $\Delta i$ from this nominal value.

Based on this information, we need to check whether $q(p1,\u2026,pn)\u2264t$ for all possible combinations of values $pi$ from the corresponding intervals $[p\u02dci\u2212\Delta i,p\u02dci+\Delta i]$.

We can perform this checking as follows:

- 1.First, we apply the algorithm $Q$ to compute the value(45)$Q\u02dc=Q(p\u02dc1,\u2026,p\u02dcn)$
- 2.Then, for each $i$ from 1 to $n$, we apply the algorithm $Q$ to compute the value(46)$Qi=Q(p\u02dc1,\u2026,p\u02dci\u22121,p\u02dci+\Delta i,p\u02dci+1,\u2026,p\u02dcn)$
- 3.Then, we compute(47)$M=Q(p\u02dc1\u2212\Delta 1,\u2026,p\u02dcn\u2212\Delta n)$
- 4.We compute(48)$Q\u02dcnew=1n+2\xb7(Q\u02dc+\u2211i=1nQi+M)$
- 5.We computewhere $k0$ depends on the level of confidence that we can achieve;(49)$b=Q\u02dcnew+\u2211i=1n|Qi\u2212Q\u02dcnew|+k0\xb72n+1\xb7\u03f53$
- 6.
Finally, we check whether $b\u2264t$.

If $b\u2264t$, this means that the desired specifications are always satisfied. If $b>t$, this means that we cannot guarantee that the specifications are always satisfied.

*Comment.*For the Cauchy method, similarly, after computing $Q\u02dc=Q(p\u02dc1,\u2026,p\u02dcn)$ and $Q(k)=Q(p\u02dc1+\eta 1(k),\u2026,p\u02dcn+\eta n(k))$, we can compute the improved estimate $Q\u02dcnew$ for $q\u02dc$ as

## Experimental Testing

### Description of the Case Study.

We tested our approach on the example of the seismic inverse problem in geophysics, where we need:

- •
To reconstruct the velocity of sound $vi$ at different spatial locations and at different depths.

- •
Based on the times $tj$ that it takes for a seismic signal to get from several setup explosions to different seismic stations.

In this example, we are interested in the velocity of sound $q$ at different depths and at different locations. To estimate the desired velocity of sound $q$, as parameters $p1,\u2026,pn$, we use travel times $tj$.

For most materials, the velocity of sound is an increasing function of density (and of strength). Thus, e.g., in geotechnical engineering, the condition that the soil is strong enough to support a road or a building is often described in terms of a requirement that the corresponding velocity of sound exceeds a certain threshold: $q\u2265t$.

*Comment.* This inequality looks somewhat different from the usual requirement $q\u2264t$. However, as we will see, the algorithm actually produces the *inverse* values $s=def1/v$. In terms of the inverse values $s$, the requirement $q\u2265t$ takes the usual form $s\u2264t0$, where $t0=def1/t$.

### Description of the Corresponding Algorithm.

As an algorithm $Q(p1,\u2026,pn)$ for estimating the velocity of sound based on the measured travel times $p1,\u2026,pn$, we used (a somewhat improved version of) the finite element technique that was originated by Hole [14]; the resulting techniques are described in Refs. [15–17].

According to Hole’s algorithm, we divide the three-dimensional volume of interest (in which we want to find the corresponding velocities) into a rectangular three-dimensional grid of $N$ small cubic cells. We assume that the velocity is constant within each cube; the value of velocity in the $i$th cube is denoted by $vi$. Each observation $j$ means that we know the time $tj$ that it took the seismic wave to go from the site of the corresponding explosion to the location of the observing sensor.

This algorithm is iterative. We start with the first-approximation model of the Earth, namely, with geology-motivated approximate values $vi(1)$. At each iteration $a$, we start with the values $vi(a)$ and produce the next approximation $vi(a+1)$ as follows.

First, based on the latest approximation $vi(a)$, we simulate how the seismic waves propagate from the explosion site to the sensor locations. In the cube that contains the explosion site, the seismic signal propagates in all directions. When the signal’s trajectory approaches the border between the two cubes $i$ and $i\u2032$, the direction of the seismic wave changes in accordance with Snell’s law $sin(\theta i)/sin(\theta i\u2032)=vi(a)/vi\u2032(a)$, where $\theta i$ is the angle between the seismic wave’s trajectory in the $i$th cube and the vector orthogonal to the plane separating the two cubes. Snell’s law enables us to find the trajectory’s direction in the next cube $i\u2032$. Once the way reaches the location of the sensor, we can estimate the travel time as $tj(a)=\u2211i\u2113ji/vi(a)$, where the sum is taken over all the cubes through which this trajectory passes, and $\u2113ji$ is the length of the part of the trajectory that lies in the $i$th cube.

### Reasonable Way to Gauge the Quality of the Resulting Estimate for the Velocity of Sound vi.

A perfect solution would be to compare our estimates with the actual velocity of sound at different depths and different locations. This is, in principle, possible: we can drill several wells that are in different locations and directly measure the velocity of sound at different depths. In practice, however, such a drilling is extremely expensive—this is why we use the seismic experiment to measure this velocity indirectly.

This is indeed a practical problem in which it is important to take model inaccuracy into account. In this problem, there are two sources of uncertainty.

The first is the uncertainty with which we can measure each travel time $tj$. The travel time is the difference between the time when the signal arrives at the sensor location and the time of the artificially set explosion. The explosion time is known with a very high accuracy, but the arrival time is not. In the ideal situation, when the only seismic signal comes from the our explosion, we could exactly pinpoint the arrival time as the time when the sensor starts detecting a signal. In real-life, there is always a background noise, so we can determine the arrival time only with some inaccuracy.

The second source of uncertainty comes from the fact that our discrete model is only an approximate description of the continuous real Earth. Experimental data show that this second type of uncertainty is important, and it cannot be safely ignored.

Thus, our case study is indeed a particular case of a problem in which it is important to take model inaccuracy into account.

### Estimating Uncertainty of the Result of Data Processing: Traditional Approach.

To compare the new method with the previously known techniques, we use the uncertainty estimate for this problem performed in Refs. [15–17], where we used the Cauchy-based techniques to estimate how the measurement uncertainty affects the results of data processing.

In accordance with this algorithm, first, we computed the value $Q\u02dc=Q(p\u02dc1,\u2026,p\u02dcn)$ by applying the above-modified Hole algorithm $Q(p1,\u2026,pn)$ to the measured travel times $q\u02dci=t\u02dci$.

After that, we simulated the Cauchy-distributed random variables $\eta i(k)$ and applied the same Hole algorithm to the perturbed values $p\u02dci+\eta i(k)$, computing the values $Q(k)=Q(p\u02dc1+\eta 1(k),\u2026,p\u02dcn+\eta n(k))$. Based on the differences $\Delta q(k)=defQ(k)\u2212Q\u02dc$, we then estimated the desired approximation error $\Delta $.

### Let us Now Apply the New Approach to the Same Problem.

In the new approach, instead of using the original value $Q\u02dc=Q(p\u02dc1,\u2026,p\u02dcn)$, we use a new estimate $Q\u02dcnew$—which is computed using Eq. (50).

Then, instead of using the original differences $\Delta q(k)=defQ(k)\u2212Q\u02dc$, we use the new differences $\Delta qnew(k)=defQ(k)\u2212Q\u02dcnew$. After that, we estimate the value $\Delta new$ by applying the maximum-likelihood method to these new differences.

### Which Estimate is More Accurate.

To check which estimates for the velocity of sound are more accurate—the estimates $Q\u02dc$ produced by the traditional method or the estimates $Q\u02dcnew$ produced by the new method—we use the above criterion for gauging the quality of different estimates. Specifically, for each of the two methods, we computed the RMS approximation error $1/n\xb7\u2211i(t\u02dci\u2212ti)2$ describing how well the travel times $tt\u02dci$ predicted based on the estimated velocities of sound match the observations $ti$.

We performed this comparison 16 times. In one case, the RMS approximation error increased (and not much, only by 7%). In all other 15 cases, the RMS approximation error decreased, and it decreased, on average, by 15%.

This result shows that the new method indeed leads to more accurate predictions.

## Future Work: Can We Further Improve the Accuracy

### How to Improve the Accuracy: A Straightforward Approach.

As we have mentioned, the inaccuracy $Q\u2260q$ is caused by the fact that we are using a finite element method with a finite size element. In the traditional finite element method, when we assume that the values of each quantity within each element are constant, this inaccuracy comes from the fact that we ignore the difference between the values of the corresponding parameters within each element. For elements of linear size $h$, this inaccuracy $\Delta x$ is proportional to $x\u2032\xb7h$, where $x\u2032$ is the spatial derivative of $x$. In other words, the inaccuracy is proportional to the linear size $h$.

A straightforward way to improve the accuracy is to decrease $h$. For example, if we reduce $h$ to $h2$, then we decrease the resulting model inaccuracy $\u03f5$ to $\u03f5/2$.

This decrease requires more computations. The number of computations is, crudely speaking, proportional to the number of nodes. Because the elements fill the original area and each element has volume $h3$, the number of such elements is proportional to $h\u22123$.

This leads to decreasing the inaccuracy by a factor of $h/h\u2032$, which is equal to $K3$.

For example, in this straightforward approach, if we want to decrease the accuracy in half $(K3=2)$, we will have to increase the number of computation steps by a factor of $K=8$.

### Alternative Approach: Description.

### Why This Approach Decreases Inaccuracy.

We know that $Q(p1+\Delta p1,\u2026,pn+\Delta pn)=q(p1+\Delta p1,\u2026,pn+\Delta pn)+\delta q$, where, in the small vicinity of the original tuple $(p1,\u2026,pn)$, the expression $q(p1+\Delta p1,\u2026,pn+\Delta pn)$ is linear, and the differences $\delta q$ are independent random variables with zero mean.

Due to linearity and the fact that $\u2211k=1K\Delta i(k)=0$, the first average in Eq. (55) is equal to $q(p1,\u2026,pn)$. The second average is the average of $K$ independent identically distributed random variables, and we have already recalled that this averaging decreases the inaccuracy by a factor of $K$.

Thus, in this alternative approach, we increase the amount of computations by a factor of $K$, and as a result, we decrease the inaccuracy by a factor of $K$.

### New Approach is Better Than the Straightforward One.

In general, $K>K3$. Thus, with the same increase in computation time, the new method provides a better improvement in accuracy than the straightforward approach.

*Comment.* The above computations refer only to the traditional finite element approach, when we approximate each quantity within each element by a constant. In many real-life situations, it is useful to approximate each quantity within each element not by a constant, but rather by a polynomial of a given order (see, e.g., [18]): by a linear function and by a quadratic function. In this case, for each element size $h$, we have smaller approximation error but larger amount of computations. It is desirable to extend the above analysis to such techniques as well.

## Acknowledgment

This work was supported in part by the National Science Foundation grants HRD-0734825 and HRD-1242122 (Cyber-ShARE Center of Excellence) and DUE-0926721. The authors are greatly thankful to the anonymous referees for valuable suggestions.