AIAA JOURNAL
Vol. 35, No. 4, April 1997
Analysis of Uncertain Structural Systems
Using Interval Analysis
S. S. Rao¤
Purdue University, West Lafayette, Indiana 47907-1288
and
L. Berke²
NASA Lewis Research Center, Cleveland,Ohio 44135
The imprecision or uncertainty present in many engineering analysis/design problems can be modeled using
probabilistic, fuzzy, or interval methods. This work considers the modeling of uncertain structural systems using
interval analysis.By representing each uncertain inputparameteras an interval number, a static structural analysis
problem can be expressed in the formof a system of linear interval equations. In addition to the direct andGaussian
elimination-based solution approaches, a combinatorial approach (based on an exhaustive combination of the
extreme values of the interval numbers) and an inequality-basedmethod are presented for ® nding the solution of
interval equations. The range or interval of the solution vector (response parameters) is found to increase with
increasing size of the problem in all of the methods. An interval± truncation approach is proposed to limit the
growth of intervals of response parameters so that realistic and accurate solutions can be obtained in the presence
of large amounts of uncertainty. Numerical examples are presented to illustrate the computational aspects of the
methods and also to indicate the importance of the truncation approach in practical problems. The utility of
interval methods in predicting the extreme values of the response parameters of structures is discussed.
I. Introduction
A LL engineering design problems involve imprecision or ap-proximationor uncertaintyto varyingdegrees.As an example,
consider a beam whose length may be stated in different ways as
follows: ª lies between 3.4 and 3.6 m,º or ª about 3.5 m,º or ª has a
mean value of 3.5 m and a standard deviation of 0.05 m, and fol-
lows normal distribution.º Depending on the nature of imprecision,
the analysis/design of the system can be conducted using interval
analysis,1 , 2 fuzzy theory,3 or a probabilistic approach.4 In interval
analysis, the uncertain parameter is denoted by a simple range. In
addition to the range, if a preferencefunction is used to describe the
desirability of using different values within the range, fuzzy theory
can be used. In fuzzy theory, the imprecision is interpreted as the
designer’s choice to use a particular value for the uncertain param-
eter. On the other hand, if the uncertain parameter is described as a
random variable following a speci® ed probability distribution, the
probabilistic approach can be used. It can be seen that when infor-
mation about an uncertain parameter in the form of a preference
or probability function is not available, the interval analysis can be
used most conveniently.
The following examples indicate typical situations in which an
uncertain parameter can be modeled as an interval number.
1) In design and manufacture, if a geometric parameter x is sub-
jected to tolerances as x § D x , then the parameter has to be treated
as an interval number as x = [x ¡ D x , x + D x].
2) In engineeringanalysis, certain parameters such as the magni-
tude of wind load (P) may be known to vary over a range P1 to P2;
but the probability distribution of P may not be known. In such a
case, the load can be treated as an interval number as P = [P1, P2].
3) In the sensitivityanalysisof a system,wemight be interestedin
® nding the in¯ uenceof changinga parameter over a speci® ed range,
x § D x . In such a case, the in¯ uence of changes in the independent
variable x on a dependent variable (response quantity) f can be
represented by an interval as f = [ f ¡ D f, f + D f ] where§D f
denotes the variations in f caused by the changes§D x in x .
Received Sept. 6, 1995;accepted for publication June18, 1996.Copyright
c° 1997 by S. S. Rao and L. Berke. Published by the American Institute ofAeronautics and Astronautics, Inc., with permission.
¤ Professor, School of Mechanical Engineering. Member AIAA.
² Structures Division. Associate Fellow AIAA.
4)Many engineeringanalysesuse a Taylor-seriesapproximation.
The absolutevalueof the error involvedin usingonly a ® nitenumber
of terms in the Taylor series is always known so that the actualvalue
of the functionbeing representedbecomesan interval. For example,
if f (x) is expanded around x ¤ using the series
Äf (x) = Äf (x ¤ + D x) ¼ f (x ¤ ) + f 0 (x ¤ ) D x + f 0 0 (x ¤ ) ( D x)
2
2!
+ ¢ ¢ ¢ + f (k) (x ¤ ) ( D x)
k
k!
(1)
the error is given by
e = f
(k + 1)( n )
( D x)k + 1
(k + 1)!
x ¤ < n < x ¤ + D x (2)
so that the actualvalue of f (x) lies in the interval f (x)= [ Äf (x) ¡ e,Äf + e].
5) The performance characteristics of most engineering sys-
tems vary during their lifetimes because of aging, creep, wear, and
changes in operatingconditions.The design of such systems should
take care of the intervalsof deviationof variablesfrom their nominal
values.
Although the analysis/design problems are to be solved using
intervals for the design parameters, most traditional methods con-
sider only nominal values for simplicity and convenience.A fuzzy-
analysis-basedapproachwas presentedby Rao and Sawyer5 for the
descriptionof systems containing information and features that are
vague, imprecise, qualitative, linguistic, or incomplete. In particu-
lar, a fuzzy ® nite element method was developed for predicting the
response of systems governedby linear systems of equations. In the
reliability analysis of structures using stress±strength interference
theory, an estimation of the failure probability of the structure is
based on a precise knowledge of the exact probability distributions
of the random parameters.6 Even then, the numerical computation
of the extreme values of the induced stress requires the modeling of
tail portionsof the probabilitydistributionsof the randomvariables,
which is very dif® cult in practice. It can be seen from the literature
that interval analysis methods were not used in engineering appli-
cations.
This paper presents methods that can be used for the analysis of
engineering systems for which the input parameters are given as
727
D
ow
nl
oa
de
d
by
N
O
RT
H
W
ES
TE
RN
P
O
LY
TE
CH
IC
A
L
U
N
IV
. o
n
A
pr
il
12
, 2
01
4
| ht
tp:
//a
rc.
aia
a.o
rg
| D
OI
: 1
0.2
514
/2.
164
728 RAO AND BERKE
interval numbers. Because most engineering systems can be ana-
lyzed using the ® nite element method, the problem considered in
this work can be stated as follows: Given a system of linear equa-
tions that describe the response of a discretized structural model in
the form
[A( p)]X = B( p) (3)
and
R = R(X) (4)
where [A] is the stiffness matrix of the system, p is the vector of
uncertain input parameters described as interval numbers as
p =
ìïïï
í
ïïïî
p1
p2
...
pm
üïïï
ý
ïïïþ
=
ìïïïï
í
ïïïïî
( p
Å
1
, Åp1)
( p
Å
2
, Åp2)
...
( p
Å m
, Åpm)
üïïïï
ý
ïïïïþ
(5)
where a bar below (above) a symbol denotes a lower (upper) bound,
X is the nodaldisplacementvectorof the discretestructuralmodel,B
is the vector of nodal forces (which might depend on the uncertain
parameters), and R is the vector of response quantities, such as
element stresses and strains, which can be expressed in terms of
the nodal displacements,X, of the structure. Because R implicitly
depends on p, there will be a range for each component of R so that
R = [R
Å
, ÅR] = f (rÅ i , Åri ) g (6)
The problem is to ® nd the range of each response quantity (r
Å i
, Åri )
due to the uncertainty present in the input parameters.
The solutionof linear interval equationsand the subsequentcom-
putationof the responseparameters,Eqs. (3) and (4), basedon direct
and Gaussian elimination-basedmethods is discussed in this work.
Alternative techniques that ® nd the bounds on the solution using a
combinatorialapproach and certain inequalitiesare proposed.It has
been observed that the widths (intervals)of the responseparameters
predictedbecomewider than the true widths with an increase in the
size of the problem.To avoid the unnecessarygrowthof the intervals
of the response parameters, a computationallyef® cient approxima-
tion procedure, termed the truncation method, is presented. The
computationalaspects and the results given by the variousmethods
are discussedwith the help of numerical examples.
II. Statement of the Problem
Most linearstructuralanalysisproblemsare solvedusingthe ® nite
elementmethodwherein the governingequationsare representedas
a system of algebraic equations:
[A]X = B (7)
where [A] = [ai j ], X = f xi g , andB = f bi g . If the systemparameters
are imprecise, the coef® cients ai j and the right-hand-sideconstants
bi in Eq. (7) are not sharply de® ned. In such a case, each of the
coef® cients and constants is represented by a closed interval as
ai j = [a
Å i j
, Åai j ]; i, j = 1, 2, . . . , n (8)
bi = [b
Å i
, Åbi ]; i = 1, 2, . . . , n (9)
Equation (7) also can be written as
([A
Å
], [ ÅA])X = (B
Å
, ÅB) (10)
The solution of Eq. (7) is understood to mean
X = {X0 j [A0]X0 = B0 ,[A0] 2 [A]= ([AÅ ], [ ÅA]), B0 2 B= (BÅ , ÅB)} (11)
where all [A0] 2 [A] are assumed to be nonsingular.Equation (11)
denotes a set of solutions (solution domain) in X space such that
for every [A0] 2 [A] and every B0 2 B, there exists an X0 that
satis® es the relation [A0]X0 = B0. Note that the solution set, given
by Eq. (11), may be empty in some cases; that is, given some [A]
and B, there may not be an X that satis® es Eq. (7).
The solution set and its hull of Eq. (7) have been investigated by
several researchers.7, 8 For systemsgovernedbyEq. (7), two typesof
problemscan be identi® ed.9 In the ® rst type of problems,all possible
vectors B are sought for speci® ed interval matrix [A] and interval
vectorX. This gives the range of B. In the second type of problems,
we seek to ® nd the interval vector X for speci® ed interval matrix
[A] and interval vector B. This gives the domain of solutionX.
The objective of this work is to investigate the solution methods
for the second type of problem. The direct approach, the Gaussian
elimination technique, a combinatorial approach, a bounding tech-
nique (based on certain inequalities),and a truncation approach are
presented for the solution of Eq. (7). Although the direct method,
the Gaussian elimination approach, and the inequalities are well
established for crisp problems, their application to interval prob-
lems, particularly to those arising in engineering analyses, is new.
Furthermore, the truncation technique presented in this work rep-
resents a completely new contribution.The numerical results given
by the variousmethods are compared for typical structural analysis
problems. Although the numerical examples are very simple, the
interval analysis methods presented are expected to be applicable
to other types of engineeringproblemswhose behavior is governed
by systems of linear equations.
III. Interval Analysis Using Direct and Gaussian
Elimination Techniques
For simple systems (where the number of equations is two or
three), the solution of Eq. (7) can be expressed in terms of the ele-
ments of the matrix [A] and the vectorB. In such cases, the response
of the system(suchas elementstresses)alsocan be determinedusing
the interval arithmetic operations given in Appendix A. However,
the widths or intervals of the responses so computed will be larger
than the actual widths.
The solution of Eq. (7) also can be found using the Gaussian
elimination method. For this, let [D] be an interval matrix whose
elements contain ranges of the corresponding elements of [A] ¡ 1 in
Eq. (7). Then the product of [D] and the interval vectorB gives the
vector X whose components represent an interval that contains the
possiblevalues of the correspondingcomponentsof the solutionsof
Eq. (7) when the elements of B are allowed to vary over the respec-
tive speci® ed ranges. When exact interval arithmetic is used in the
computations,we obtain an intervalmatrix, [A]¡ 1, whose elements
are wider than the correspondingranges of the exact inversematrix.
The computations are shown for a sample problem in Appendix B.
There is another practical dif® culty associated with the imple-
mentation of the Gaussian elimination method for the solution of
interval equations. This is related to reducing the coef® cients of
off-diagonal variables to zero during pivot operations. The basic
arithmetic operations, related to interval numbers and interval ma-
trices, are summarized in Appendix A. Because the subtraction of
an intervalnumber from itself does not result in zero (for example, if
A = [a
Å
, Åa], then A ¡ A = [aÅ , Åa] ¡ [aÅ , Åa] = [aÅ ¡ Åa, Åa ¡ aÅ ] 6= [0, 0]),the Gaussian elimination cannot be implemented accurately.
IV. Analytical Developments
A. CombinatorialApproach
If f (x1 , x2 , . . . , xn ) denotes an arbitrary (monotonic) functionof
the impreciseparametersx1, x2 , . . . , xn , the rangeof f can be deter-
mined by consideringall possible combinationsof the endpoints of
the imprecise parameters. For this, let the parameter xi be denoted
as an interval number as
xi = [x
(1)
i , x
(2)
i ] = [xÅ i
, Åxi ]; i = 1, 2, . . . , n (12)
Then all possible values of f can be determined as
fr = f (x
(i)
1 , x
( j )
2 , . . . , x
(k)
n
)
i = 1, 2; j = 1, 2; ¢ ¢ ¢ ; k = 1, 2; r = 1, 2, . . . , 2n (13)
where fr denotes the value of f for a particular combination of the
endpoints of the intervals of x1, x2, . . . , xn . The function f can be
represented as an interval number as
f = [ f
Å
, Åf ] = [ min
r = 1,2,..., 2n
( fr ), max
r = 1,2,..., 2n
( fr )] (14)
D
ow
nl
oa
de
d
by
N
O
RT
H
W
ES
TE
RN
P
O
LY
TE
CH
IC
A
L
U
N
IV
. o
n
A
pr
il
12
, 2
01
4
| ht
tp:
//a
rc.
aia
a.o
rg
| D
OI
: 1
0.2
514
/2.
164
RAO AND BERKE 729
Equations (3) and (4) are based on linear structural analysis and
hence Eq. (14) always gives the correct extreme values of the
response parameters such as nodal displacements and element
stresses.
B. Bounds on Solution Using Inequalities
Any system of interval equations of the form
[A]X = B (15)
with [A] = (([A]
Å
, [ ÅA])) and B = ((B
Å
, ÅB)) can be rewritten as
(([A0] ¡ [D A], [A0]+ [D A]))X = ((B0 ¡ D B, B0 + D B)) (16)
where
[A0] = 12 (([AÅ
] + [ ÅA])), [ D A] = 12 (([ ÅA] ¡ [AÅ ])) ¸ [0] (17)
B0 = 12 (BÅ +
ÅB), D B = 12 ( ÅB ¡ BÅ ) ¸ 0 (18)
It can be observed that 1) ¡ [ D A]j X j is always less than or equal to
both ¡ [D A]X and [ D A]X and 2) [D A] j X j is always greater than or
equal to both ¡ [D A]X and [D A]X, where
j X j =
ìïïï
í
ïïïî
j x1 j
j x2 j
...
j xn j
üïïï
ý
ïïïþ
The satisfaction of the following inequality is necessary and suf® -
cient for X to be a solution of Eq. (15) (Refs. 10 and 11):
[D A]j X j + D B ¸ j [A]X ¡ B j (19)
The inequality (19) is equivalent to
[A]X ¡ [D A] j X j · ÅB (20)
and
[A]X + [D A] j X j ¸ BÅ (21)
All vectors X that satisfy the inequalities (20) and (21) constitute
the solution domain of Eq. (15). Equations (20) and (21) denote a
set of 2n inequalities,which can be expressed as
D aTi j X j + D Bi ¸ aTi X ¡ Bi ; i = 1, 2, . . . , n (22)
and
D aTi j X j + D Bi ¸ ¡ (aTi X ¡ Bi ) ; i = 1, 2, . . . , n (23)
or
l+i (X ) = a
T
i X ¡ Bi ¡ D aTi j X j ¡ D Bi ·0; i = 1, 2, . . . , n
(24)
and
l ¡i (X ) = aTi X ¡ Bi + D aTi j X j + D Bi ¸ 0; i = 1, 2, . . . , n
(25)
where Bi and D Bi denote the i th components ofB and D B, respec-
tively, ai and D ai represent the i th rows of the matrices [A] and
[ D A], respectively,and the superscriptT indicates the transpose. It
can be seen that the vectorX is restricted to the convex polyhedron
bounded by the hyperplanes l+i = 0 and l ¡i = 0. Because there are
only n unknowns in the vector X, at any vertex of the polyhedron,
the equality sign holds only for n of the 2n conditions of Eqs. (24)
and (25). It can be observed that the hyperplanes l+i = 0 and l ¡i = 0
lie symmetrically with respect to the hyperplaneaTi X ¡ Bi = 0. As
such, they cannot intersect in the orthant containing X0 where X0
satis® es the equations [A0]X0 = B0 (Ref. 10). This implies that
l+i = 0 and l ¡i > 0 (26)
or
l+i < 0 and l ¡i = 0 (27)
Thus the 2n vertices of the polyhedronare characterizedby the fact
that the equality sign holds only for one of the equations (24) and
(25).
If all of the admissible vectors X do not lie in the same orthant,
the solution of Eq. (15) is given by the union of convex polyhedra,
which may not be convex. To determine the vertices of the solution
set of Eq. (15), we need to ® nd the minimum/maximum values of
the various components x j of the solution. This can be achieved by
solving the following linear programming problems10:
minimize/maximize x j
subject to
l+i (X) ·0; i = 1, 2, . . . , n
¡ l ¡i (X) ·0; i = 1, 2, . . . , n
(28)
It can be seen that the problemof Eq. (28) with minimization (max-
imization) gives the lower (upper) boundon the interval of x j . Thus
the solution of 2n linear programming problems provides the solu-
tion of Eq. (15).
C. Interval Analysis with Truncation
From the de® nitions of interval arithmetic operations, given by
Eq. (A1), it can be seen that the evaluationof any rationalexpression
in interval arithmetic will lead to the exact range of values of the
real rational function de® ned by the expression if each variable
occurs only once in the expression.1 Because it is not possible to
rewrite complicated and/or implicit expressions in which a number
of interval variables appear several times, there is no simple way to
obtain the exact intervals of functions involved in most engineering
applications. Thus it can be expected that the width of the interval
of a function grows with the number of interval variables and the
number of arithmetic operations (or size of the problem) involved.
To limit the growth of intervals of response parameters for large
amounts of uncertainties, a truncation procedure is proposed. This
procedurecan be describedas follows.Let a = [a
Å
, Åa] and b = [b
Å
, Åb]
denote the input interval variables and c = [c
Å
, Åc] the result of an
interval arithmetic operation. If the central values of the variables
a0 = (a
Å
+ Åa)/ 2 and b0 = (b
Å
+ Åb)/ 2 are used, the corresponding
arithmetic operation on the central numbers leads to the result c0.
Let c0 be called the crisp (or central) result that can be obtained
when the variability is absent in the input variables. If c0 is close to
zero, no truncation is used. On the other hand, if c0 is not close to
zero, then the relative deviation ( D ) of the interval, Åc ¡ cÅ , from thecentral value is computed as
D = D 1 + D 2 (29)
where
D 1 =
ê
ê
ê
ê
c
Å ¡ c0
c0
ê
ê
ê
ê
, D 2 =
ê
ê
ê
ê
c
Å
0 ¡ Åc
c0
ê
ê
ê
ê
(30)
Because the deviation D is expected to be larger than the true devi-
ation, the maximum permissible relative width of c is speci® ed as
2t . The value of t can be speci® ed easily on the basis of the known
deviations of the input variables a and b. For example, the value of
t can be speci® ed to be equal to the maximum of the relative devi-
ations of the input parameters a and b from their respective central
values. The interval number c = [c
Å
, Åc] is truncated to obtain the
result c ¼ [dÅ , Åd] as follows:
d
Å
= c
Å
, Åd = Åc if D ·2t (31)
(if both D 1 and D 2 are smaller than or equal to the permissible
deviation t, no truncation is used),
d
Å
= c0 + t (c
Å ¡ c0), Åd = c0 + t ( Åc ¡ c0) if D > 2t (32)
(if both D 1 and D 2 are larger than the permissible deviation t,
[(c ¡ c0)/c0] is truncated at ¡ t and +t),
d
Å
= c
Å
, Åd = c0 + t ( Åc ¡ c0) if D 1 · t and D 2 > t (33)
D
ow
nl
oa
de
d
by
N
O
RT
H
W
ES
TE
RN
P
O
LY
TE
CH
IC
A
L
U
N
IV
. o
n
A
pr
il
12
, 2
01
4
| ht
tp:
//a
rc.
aia
a.o
rg
| D
OI
: 1
0.2
514
/2.
164
730 RAO AND BERKE
D · 2t
D > 2t
D 1 · t and D 2 > t
D 1 > t and D 2 · t
Fig. 1 Truncation procedure.
(if D 1