COMPUTATIONAL ISSUES IN FORM WITH UNIFORM RANDOM VARIABLES

. The First Order Reliability Method is well accepted as an efficient way to solving structural reliability problems with linear or moderately non-linear limit state functions. High non-linearity is introduced in reliability problems by non-linear mechanical responses, but also by correlation between the random variables and by highly non-Gaussian probability distribution functions. Correlation and non-Gaussian distributions introduce non-linearities in the mapping to the standard Gaussian space, hence making search for the design point more challenging. This article discusses computational issues in finding the design point in problems involving uniform and other bounded random variables. The discussion covers the Principle of Normal Tail Approximation, and the mapping to standard Gaussian space. It addresses three different techniques to impose the domain of bounded random variables, when mapping them back from standard Gaussian to original design space. A simple but challenging academic problem is presented, involving a simply-supported beam subject to a concentrated load of random intensity. The concentrated load can occupy a random position over the beam, following an uniform distribution. Although the underlying mechanical problem is very simple, the uniform distribution introduces severe non-linearities, which makes finding the design point a very demanding task. The article also addresses algorithms that can be used to find the design point in highly non-linear problems. It is shown that the Hassofer-Lind-Rackwitz-Fiessler (HLRF) algorithm fails to converge. The performance of the improved HRFL (iHLRF) and of the Sequential Quadratic Programming algorithms are investigated with respect to two alternative mappings to standard Gaussian space and to three techniques to impose the bounds in the inverse probability mapping.


INTRODUCTION
The First Order Reliability Method is well accepted as an efficient way to solve structural reliability problems with linear or moderately non-linear limit state functions.High nonlinearity is introduced in reliability problems by non-linear mechanical responses [6], but also by correlation between the random variables and by highly non-Gaussian probability distributions.Correlation, non-Gaussianity and generally bounded distributions introduce nonlinearities in the mapping to standard Gaussian space, making the search for the design point more demanding.
This article discusses challenges in finding the design point in problems involving uniform and other bounded random variables.The discussion covers the Principle of Normal Tail Approximation [2,3], and the mapping to the standard Gaussian space.It is observed that the conventional mapping becomes highly non-linear as, during interactive design point search, some variable is pushed against its lower or upper limit.An alternative mapping is proposed, in order to make computations more stable.As standard Gaussian space is unbounded, during design point search any bounded random variable could try to exceed its limits, causing breakdown of the computations.Hence, bounds have to be imposed when mapping random variables back to the original design space.Three different schemes are presented in this article to impose the limits of bounded random variables: bisection, reflection and limit.The three schemes are tested in combination with the traditional and the alternative (proposed) probability mapping.A challenging academic problem is also presented in the article.It involves a simply-supported beam subject to a concentrated load of random intensity.The load can occupy a random position over the beam, following an uniform distribution.Although the underlying mechanical problem is very simple, the uniform distribution introduces strong non-linearities, which makes finding the design point a challenging task.
It is shown that the Hasofer-Lind-Rackwitz-Fiessler (HLRF) algorithm fails to converge for most configurations of the proposed problems.Hence, the improved HRFL (iHLRF) and Sequential Quadratic Programming (SQP) algorithms are also employed in the solutions.The performance of these algorithms is investigated with respect to two mappings to standard Gaussian space and to the three schemes to impose the bounds in the reverse probability mapping.In total, twenty-four configurations of the proposed benchmark problem are investigated.The article is finished with relevant conclusions with respect to the several alternatives proposed herein for the solution of reliability problems involving uniform and other bounded random variables.

The First Order Reliability Method (FORM)
Let X be a random variable vector describing the uncertainties in geometry, material properties and loading.A limit state function X) is written in such a way as to divide the random variable domain in safety and failure domains: The failure probability is given by: where denotes the joint probability density of random vector .The First Order Reliability Method (FORM) is generally accepted as an efficient way to solving Eq. ( 2), for low dimensions of vector X and for linear or moderately non-linear limit state functions X).In the FORM method, Eq. ( 2) is solved indirectly, by introducing a convenient mapping from the space of the original random vector ( -space) to the so-called standard Gaussian space ( -space) [3,9,10]: where all components of vector are independent and identically distributed standard Gaussian random variables.This mapping can be accomplished by means of the Principle of Normal Tail Approximation (to be described) and by the Nataf transformation.
In standard space, the joint probability density is rotationally symmetric: hence, the point * on the limit state function which is closest to the origin represents the most probable failure point, also known as design point.This feature allows search for the design point to be cast as a constrained optimization problem: * arg min • subject to 0 From Eq. ( 4), ‖ * ‖ * • * / is the so-called Hasofer-Lind reliability index [5,11] which comes to be the distance between * and the origin of standard Gaussian space.The FORM method, therefore, consists in finding the design point * and approximating the original limit state function by a tangent hyper-surface at the design point.Hence, the firstorder approximation of the failure probability becomes: where Φ . is the cumulative standard Gaussian probability distribution function.
For the highly non-linear problems to be considered herein, the first order approximation of the failure probability is unlikely to be accurate.However, it is still important to be able to find the design point, for instance, in order to perform Monte Carlo simulation with importance sampling using the design point [4] or to construct precise response surface approximations centered on the design points [7].
It is important to point out that the concept of a design point as the most probable failure point is theoretical, and associated to the transformation to standard Gaussian space.For instance, for problems involving only uniform random variables, all points in the failure domain are equally likely in -space… For problems involving a single uniform random variable, multiple equally likely failure points would exist in -space.However, the concept and the usual interpretations remain valid in -space.

Traditional normal tail approximation
The Principle of Normal Tail Approximation [3] allows a non-Gaussian variable to be transformed in an "equivalent" Gaussian variable by means of a point-wise probability mapping: * * Eq. ( 6) is a probability-preserving equation, which makes the probability content to the "left" and "right" of point * to remain unaltered in the mapping.However, two parameters ( and -mean and standard deviation) of the "equivalent" Gaussian variable have to be determined.Hence, an additional equation is required.The second equation that has been traditionally used to complete the mapping is: * * This equation is completely arbitrary, and is not required for probability preservation.Other equations and rules for the probability mapping have also been presented elsewhere (Chen and Lind, 1983).
Introducing the Hasofer-Lind transformation [5]: the mapping can be accomplished directly to standard Gaussian space.Manipulating Eqs. ( 6) to (8), well-known equations are obtained for the desired parameters ( and ) : * Φ * * * * * Now, for illustration purposes, assume that random variable follows an uniform distribution with parameters : Figure 1 illustrates how the traditional mapping behaves as the point * approaches the upper limit ( * with → 0).It can be observed that → 0, which could cause breakdown of the numerical computations.The observed behavior also introduces strong non-linearities in the probability mapping, which can make the design point search unstable, if any of the * points comes too close to a random variable limit.approaches the upper limit of a bounded random variable ( * with → 0).

Alternative normal tail approximation
In an attempt to avoid potential problems with the traditional probability mapping, as observed in Figure 1, an alternative mapping is proposed herein.As Eq. ( 7) is arbitrary, it is proposed to use the simple rule instead: The standard deviation is obtained directly from Eq. ( 8), using the relation * Φ * , which originates in the probability-preserving equation (Eq.6). Figure 1 shows that the resulting mapping is more stable as a point * approaches the upper limit ( * with → 0).When point * is equal to the mean or mode, the alternative mapping cannot be used, and the traditional mapping is used instead (only for one interation).The alternative mapping proposed herein is compared with the traditional mapping in connection to some example problems in the sequence.

Imposing random variable bounds in the inverse mapping
With the mapping to standard Gaussian space, all variables become unbounded.Hence, when a new, trial point is evaluated during iterative design point search, in -space, it could become unfeasible for exceeding a limit of the original random variable in -space.Hence, when mapping a new point back to original design space, the random variable limits have to be imposed.Although this appears to be an obvious problem with the mapping to standard Gaussian space, it has not been explicitly addressed in the published literature.
In this section, three alternative schemes are presented for bounding limited random variables: bisection, reflection and limit.Assume that the lower and upper bounds of random variable are and , respectively.The bisection bounding scheme is: The reflection scheme is obtained by: As stated in Eq. ( 13), the reflection scheme is not guaranteed to generate feasible points.Hence, after performing the reflection, the point should be checked again using the limit scheme.
For the limit scheme, one should choose a limiting value which is compatible with the accuracy of the numerical representation of the standard Gaussian cumulative probability distribution function, Φ ., and its inverse Φ .-see first line of Eq. ( 9).In this article, 10 is considered.This corresponds to a reliability index of Φ 10 5.6.This is still safely below the operational limit of the StRAnD software [1], beyond which the numerics used in the computations crumble ( 8, Φ 8 10 ).The limit bounding scheme is simply : The three bounding schemes allow the computation to continue safely, and also allow the sequence of points to return naturally to the interior of the feasible domain, if the design point is located there.

The Hasofer, Lind, Rackwitz and Fiessler (HLRF) algorithm
The most popular algorithm used to solve the reliability problem stated in Eq. ( 4) is the Hasofer and Lind [5], Rackwitz and Fiessler [11] algorithm (HLRF).This algorithm is known to be very efficient when it converges, but is also known to fail to converge very often [8,13].The algorithm is based on finding the zero of successive linearizations of the limit state surface, leading to the well-known sequence: In Eq. ( 15), is the gradient vector of the limit-state function with respect to evaluated at .

The improved Hasofer, Lind, Rackwitz and Fiessler (iHLRF) algorithm
An improved Hasofer, Lind, Rackwitz and Fiessler (iHLRF) algorithm was presented by Zhang and Kiureghian [13].In the iHLRF algorithm, a linear search is performed: in the direction given by the HLRF algorithm: In Eqs. ( 16) and ( 17), is a search direction and α is a step size, selected from minimization of an appropriate merit function.The following merit function is proposed in ref. [13]: where is a penalty parameter.The linear search is defined as: The linear search can be performed approximately using the popular Armijo rule [13].An appropriate step α is selected, such as to obtain a sufficient decrease in the merit function: Following Zhang and Kiureghian [13], direction given in Eq. ( 17) is a descent direction of the merit function in Eq. ( 18), for ∀ ∈ , as long as the following condition is satisfied: This condition is used as a rule for updating the penalty parameter: Literature-suggested values [8,13] for the algorithmic constants are 2, 0.1 and 0.5.The value 2 is used throughout the article in the numerical examples.The recommend values 0.1 and 0.5 are also used for the Armijo line search in the numerical examples.However, is was found that the iHLRF does not converge for many of the tested problems using the recomended values.Hence, the iHLRF algorithm is also tested using 0.5 and 0.5.For convenience, the two resulting algorithms are refered to as iHLRF15 and iHLRF55, respectively, in the results section.

Sequential Quadratic Programming (SQP)
The constrained optimization problem stated in Eq. ( 4) can be converted in an unconstrained problem by introducing the associated Lagrangian problem [8,12]: If * is a solution to Eq. ( 4), then there exists a Lagrangian multiplier λ * such that [8,12]: This can be used to cast the optimization problem in terms of the following set of simultaneous equations: The Newton method can be used to solve the above system of non-linear equations iteratively: The solution to this problem is given by : This result happens to be the solution and the Lagrangian multiplier for the following quadratic optimization problem: * , λ * arg min 1 2 , λ subject to (28) Therefore, this algorithm is called Sequential Quadratic Programming algorithm (SQP).The present implementation of SQP follows the guidelines set up by Schittkowski [12], as detailed in [8].The step in Eq. ( 26) is computed using the augmented Lagrangian function : where c is a penalty parameter, also updated iteratively.
The solution in Eq. ( 27) requires the Hessian of the Lagrangian function ( , λ ), which is generally not known.In this implementation, the Hessian is replaced by an approximate matrix, which is updated iteratively, following the BFGS scheme.However, due to our annoiance with non-convergence of SQP for some problems, the exact Hessian was also used.The computations using the exact Hessian are identified as SQPeH in the numerical results session.

A SIMPLE BUT CHALLENGING BENCHMARK PROBLEM
Consider a simply-supported beam subject to a concentrated load, as illustrated in Figure 2. The load of intensity P can occupy any position over the beam spam L. Both load position and intensity are assumed to follow uniform distributions: The load position limits are 0 for all problems, and varies following Table 1.The different values of change the problem considerably, as can be observed in the design points and reliability index results presented in Table 1.Three different cases of correlation are also considered: no correlation, positive (ρ=0.5) and negative (ρ=-0.5)correlation between variables and .Also, variants of the two-dimensional problem are obtained for the load intensity following Normal and Log-normal distributions, following Table 1.In total, 24 two-dimensional problems (variants) are considered.The limit state function for the two-dimensional problems is: where the resistance is given by 4 ⁄ . Parameter controls the size of the failure domain and hence the magnitude of the failure probability.For most (conventional) problems, =0.05 is considered (Table 1).A problem similar to the one-dimensional problem above (called "corner" problem herein), is obtained for 10 .For the corner problems, the convergence criterion is not used.Variables , , , , and are treated as nondimensional.For all problems, 1.

RESULTS
Results are presented in this section for a combination of four optimization algorithms (HLRF, iHLRF with two sets of parameters, SQP), three schemes for imposing random variable bounds (bisection, reflection and limit) and for two probability mappings (traditional and alternative normal tail approximation).The combination results in 24 algorithms, which are employed in the solution of 24 two-dimensional problems and four In the problem solutions, it is observed that the alternative mapping leads to similar trajectories for close-by initial points, whereas the traditional mapping leads to more irregular trajectories.Difficulties were observed in reaching the prescribed tolerance (‖ ‖ 10 ), as towards the end of the convergence history several (successive) points remain clogged.
Figure 3 illustrates the convergence trajectories for one initial point but different optimization algorithms, also for problem 1 (Table 1).In this figure, it can be observed that the HLRF algorithm gets trapped in a circle, bouncing from one side of the mid-spam to the other, hence failing to converge.This happens for most of the 24 2D problems, when the HLRF algorithm is used.The iHLRF15, iHLRF55 and SQP algorithms do converge for this problem, as shown.The performance of different algorithms can be evaluated from Table 2, which summarizes results for the 12 2D problems with uniform probability distributions for load .As the particular starting points for which convergence was achieved are of no general interest, results are summarized for each algorithm and mapping.Details of starting point and bounding schemes are omitted.In Table 2, a number indicates that convergence was achieved for the six starting points and for the three bounding schemes (18 runs); in this case the number is the average number of limit state function calls required for convergence.Keywords are as indicated in the footnote to Table 2.
For the problems with 1.0, both positive and negative correlations lead to greater convergence difficulties.For the other problems, correlations actually speeded up convergence, leading to smaller number of limit state function calls.
In the solution of the 12 problems summarized in Table 2, the iHLRF55 was the most efficient algorithm.It converged for all problems using the traditional mapping, and for most problems using the alternative mapping.The iHLRF15 algorithm converged for fewer problems, and when it converged, it took on average 13% more limit state function calls then the iHLRF55 algorithm.The SQP algorithm also converged for all problems using the traditional mapping, but in general required from 30 to 42% more limit state function calls than the iHLRF55 algorithm.
In Table 2 it can also be observed that the proposed alternative mapping leads to convergence difficulties for the so-called corner problems.For all other problems, however, the alternative mapping leaded to an average reduction of around 10% in the required number of limit state function calls.Table 2. Summary of results for 2D Problems 1 to 12.
HLRF Keywords none and some are obvious; few means that convergence was obtained for up to three starting points (for one or more bounding schemes); most means that convergence was obtained for all but up to three starting points; mean means that convergence was obtained only for the mean starting point; an * indicates that convergence was only obtained for the limit bounding scheme.→ 1.0 0.6 0.5 1.0 0.6 0.5 1.0 0.6 0.5 1.0 0.6 0.5 Alter.
HLRF For 2D Problems with 1.0 or 0.6, the convergence sequences never tried to leave the random variable limits, for the four tested algorithms.Hence, bounding schemes were not activated.For the conventional problem with 0.5 (problems 3, 6 and 9), bounding schemes were requested occasionally, for some of the starting points.
As expected, for the corner problems the bounding schemes were more active, as the convergence sequence leaded to more frequent boundary hits.Figure 4 illustrates convergence difficulties using the alternative mapping and the limit bounding scheme.For the HLRF and iHLRF15 algorithms, the sequence remains trapped between the limits and .For the traditional mapping (Figure 5), the convergence is much smoother: the sequence "sticks" at the upper limit and follows on until the "design point" is found.For the traditional mapping, even the HLRF algorithm converges to the design point of the corner problem.This failure of the alternative mapping in allowing the design point of corner problems to be found is a little surprising.However, the sucess of the traditional mapping is certainly a consequence of the robust computational scheme implemented herein.
Convergence difficulties were encountered using the bissection and reflection bounding schemes, even when the traditional mapping is used.For the HLRF and iHLRF15 algorithms, the sequences also remain trapped near the design point, but are not able to get there.Only the SQP algorithm was able to find the design point using the bissection and reflection schemes.For all other algorithms, only the limit scheme allowed the design point to be found.Clearly, this is a particularity of the corner problems considered here, as the design point was arbitrarily located at the upper limits of the two random variables.This would not occur very often in practical problems.In general, it was also observed that the alternative mapping leads to more frequent boundary hits, in comparison to the traditional mapping.This justifies the observed convergence difficulties for the corner problem and the alternative mapping.Moreover, both iHLRF algorithms tested herein lead to more frequent boundary hits than the SQP algorithm, with appears to lead to more well-behaved convergence sequences.
Results obtained for problems 12 to 24, with Normal or Log-normal distribution of load intensity are presented in Table 3.None of the studied algorithms was able to identify the design point for all the starting points of all problems.The SQP algorithm failed to converge for most problems involving small coefficients of variation of the load (c.o. v. 0.05 for problems 13 to 15 and 19 to 21). Figure 6 illustrates convergence difficulties for SQP algorithm and problem 13.The SQP algorithm remains trapped in a circle, bouncing from one side of the mid-spam to the other.This behavior also happened for the computation using the exact Hessian (SQPeH).
The iHLRF15 and iHLRF55 algorithms solved most problems, but failed to converge for some of them.For those problems were both algorithms converged, the iHLRF55 algorithm required between 20 to 30% less limit state function calls than the iHLRF15 algorithm.

CONCLUSIONS
This article presented and discussed strategies for finding the design point in problems involving uniform and other bounded random variables.The article introduced a challenging benchmark problem, involving one concentrated load acting on a simply supported beam.Although mechanically simple, finding the design point was shown to be quite challenging, due to the uniform distribution of load positions and due to random load intensities.Several topics were discussed and illustrated with respect to results obtained for this benchmark problem.
In particular, the article discussed alternatives to be considered when random variables approach their upper or lower limits during design point search.An alternative transformation was suggested, in order to smoothen the mapping to standard Gaussian space: the mean of the equivalent normal distribution is simply given by the mean, or the mode, of the original distribution.This alternative procedure has indeed smoothened some of the studied problems, leading to convergence for some cases where the traditional normal tail approximation failed, but also reducing in around 10% the required number of limit state function calls.However, the suggested alternative mapping failed when the design point was located very close to a random variable limit.
Three schemes were presented herein to impose the bounds of the random variables in the inverse probability mapping: bisection, reflection and limit.The bisection and reflection schemes resulted in similar behavior for most problems, with the reflection scheme leading to convergence failure in a few cases.The limit scheme was shown to perform better in the (unlikely) case of a design point lying very close to a random variable limit.Indeed, for some of the presented "corner" problems, the limit scheme was the only scheme leading to convergence.
Four algorithms were studied with respect to their robustness and efficiency in finding the design point, in combination with the alternative and traditional normal tail approximation, and with the three bounding schemes tested herein.The HLRF algorithm converged for only a few cases, remaining trapped in some sort of circle for most of the problems studied.The Sequential Quadratic Programming (SQP) algorithm converged for most of the problems studied, but failed to converge for the 2D problems involving Normal and Log-normal load intensities and small c.o.v..The improved HLRF (iHLRF) algorithm converged for some of the problems studied, using the Armijo line search with the literature-recommended values of a=0.1, b=0.5.However, using the parameters a=0.5, b=0.5, the algorithm converged for a larger number of problems with a significant reduction in the mean number of limit state function calls.
Results were also obtained for multi-dimensional versions of the 2D problems presented herein.These results will be presented at the conference.

Figure 1 .
Figure 1.Traditional (left) and alternative (right) standard Gaussian mapping as point * approaches the upper limit of a bounded random variable ( * with → 0).

Figure 2 .
Figure 2. Simply-supported beam subject to a randomly-located concentrated load of random intensity.

Figure 3 .
Figure 3. Converge trajectories for different algorithms with traditional mapping; Problem 1.

Figure 4 .
Figure 4. Converge trajectories for different algorithms with alternative mapping and limit bounding; 2D Problem 12.

Figure 5 .
Figure 5. Converge trajectories for different algorithms with traditional mapping and limit bounding; 2D Problem 12.

Table 1 .
Variable data and design point results for two-dimensional problems.
8 and 16-dimensional problems.Due to strong problem non-linearities, several initial points have to be considered for each problem and algorithm.Six different starting points are considered for each problem, totaling 6048 problem solutions.Initial points considered for the two-dimensional problems

Table 3 .
Summary of results for 2D Problems 12 to 24.