Fuzzy parametric uncertainty analysis of linear dynamical systems A surrogate modeling approach

13 Pages·2012·0.43 MB·English

by  R. Chowdhury

MechanicalSystemsandSignalProcessing32(2012)5–17
ContentslistsavailableatSciVerseScienceDirect
Mechanical Systems and Signal Processing
journal homepage: www.elsevier.com/locate/ymssp
Fuzzy parametric uncertainty analysis of linear dynamical systems:
A surrogate modeling approach
R. Chowdhury, S. Adhikarin
CollegeofEngineering,SwanseaUniversity,SingletonPark,SwanseaSA28PP,UK
a r t i c l e i n f o a b s t r a c t
Articlehistory: Uncertainty propagation engineering systems possess significantcomputational chal-
Received13May2011 lenges.Thispaperexploresthepossibilityofusingcorrelatedfunctionexpansionbased
Receivedinrevisedform metamodellingapproachwhenuncertainsystemparametersaremodeledusingFuzzy
27April2012 variables. In particular, the application of High-Dimensional Model Representation
Accepted3May2012
(HDMR)isproposedforfuzzyfiniteelementanalysisofdynamicalsystems.TheHDMR
Availableonline2June2012
expansion is a set of quantitative model assessment and analysis tools for capturing
Keywords: high-dimensionalinput–outputsystembehaviorbasedonahierarchyoffunctionsof
Fuzzysystems increasing dimensions. The input variables may be either finite-dimensional (i.e., a
Structuraldynamics vector of parameters chosen from the Euclidean space RM) or may be infinite-
Uncertaintypropagation dimensional as in the function space CM½0;1(cid:2). The computational effort to determine
HDMR
the expansion functions using the alpha cut method scales polynomially with the
numberofvariablesratherthanexponentially.Thislogicisbasedonthefundamental
assumption underlying the HDMR representation that only low-order correlations
amongtheinputvariablesarelikelytohavesignificantimpactsupontheoutputsfor
most high-dimensional complex systems. The proposed method is integrated with a
commercialFiniteElementsoftware.Modalanalysisofasimplifiedaircraftwingwith
Fuzzyparametershasbeenusedtoillustratethegeneralityoftheproposedapproach.In
the numerical examples, triangular membership functions have been used and the
resultshavebeenvalidatedagainstdirectMonteCarlosimulations.
&2012ElsevierLtd.Allrightsreserved.
1. Introduction
Theconsiderationofuncertaintiesinnumericalmodelstoobtainthevariabilityofvibrationresponseisbecomingmore
desirableforindustrialscalefiniteelementmodels.Broadlyspeaking,therearetwoaspectstothisproblem.Thefirstisthe
quantification of parametric and/or non-parametric uncertainties associated with the model and the second is the
propagation of uncertainties through the model. This paper is concerned with the second problem. Several types of
uncertainties can exist in a physics based computational model. Such uncertainties include, but are not limited to (a)
parameteruncertainty(e.g.,uncertaintyingeometricparameters,frictioncoefficient,strengthofthematerialsinvolved);
(b)modeluncertainty(arisingfromthelackofscientificknowledgeaboutthemodelwhichisaprioriunknown);and(c)
experimental error (uncertain and unknown errors percolate into the model when they are calibrated against
experimentalresults).Theseuncertaintiesmustbeassessedandmanagedforcrediblecomputationalpredictions.
nCorrespondingauthor.Tel.:þ441792602088;fax:þ441792295676.
E-mailaddresses:[email protected](R.Chowdhury),[email protected](S.Adhikari).
0888-3270/$-seefrontmatter&2012ElsevierLtd.Allrightsreserved.
http://dx.doi.org/10.1016/j.ymssp.2012.05.002
6 R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17
When substantial statistical information exists, the theory of probability and stochastic processes offer a rich
mathematical framework to represent such uncertainties. In a probabilistic setting, the data (parameter) uncertainty
associatedwiththesystemparameters,suchasthegeometricpropertiesandconstitutiverelations(i.e.,Young’smodulus,
massdensity,Poisson’sratio,dampingcoefficients),canbemodeledasrandomvariablesorstochasticprocessesusingthe
so-called parametric approach. These uncertainties can be quantified and propagated, for example, using the stochastic
finite element method [1–4]. The reliable application of probabilistic approaches requires information to construct the
probability density functions of uncertain parameters. This information may not be easily available for many complex
practicalproblems.Insuchsituations,non-probabilisticapproaches,suchasintervalalgebra[5],convexmodels[6]and
fuzzy set [7] based methods, can be used. In this paper the uncertain variables describing the system parameters are
modeledusingfuzzyvariables.
Overthepastdecade,significantdevelopmentshavetakenplaceinfuzzyfiniteelementanalysisoflineardynamical
systems (see, for example, the review papers [8–10]). A fuzzy variable can be viewed as a generalization of an interval
variable.Whenanuncertainvariableismodeledusingtheintervalapproach,thevaluesofthevariableliewithinalower
andanupperbound.Thefuzzyapproachgeneralizesthisconceptbyintroducingamembershipfunction.Inthecontextof
computational mechanics, the aim of a fuzzy finite element analysis is to obtain the range of certain output quantities
(suchasdisplacement,accelerationandstress)giventhemembershipofdatainthesetofinputvariables.Thisproblem,
knownastheuncertaintypropagationproblem,hastakencenterstageinrecentresearchactivitiesinthefield.Inprinciple,
anuncertaintypropagationproblemcanbealwayssolvedusingtheso-calleddirectMonteCarlo simulation.Usingthis
approach,alargenumberofmembersoftheparametersetareindividuallysimulatedandboundsareobtainedfromthe
resultingoutputs.Formostpracticalproblems,directMonteCarlosimulationisprohibitivelycomputationallyexpensive.
Therefore, the aim of the majority of current research is to reduce the computational cost. Since fuzzy variables are a
generalization of interval variables, methods applicable for interval variables, such as classical interval arithmetic [5],
affine analysis [11,12] or vertex theorems [13], can be used. The Neumann expansion [14], the transformation method
[15–18], and more recently, response surface based methods [19,20] have been proposed for fuzzy uncertainty
propagation.Inthecontextofdynamicalsystems,severalauthorshaveextendedclassicalmodalanalysistofuzzymodal
analysis[21–27].
Acrucialfactorintheefficiencyofallthemethodsisthenumberofinputvariables,oftenknownasthe‘dimensionality’
of the problem. Here we investigate the use of the High-Dimensional Model Representation (HDMR) approach [28] to
addressthisissue.TheHDMRexpansionisanefficientformulationforhigh-dimensionalmappingincomplexsystemsif
thehigherordervariablecorrelationsareweak,therebypermittingtheinput–outputrelationshipbehaviortobecaptured
by the terms of low-order. The computational effort to determine the expansion functions using the a-cut method will
scale polynomically with the number of variables rather than exponentially. This logic is based on the fundamental
assumptionunderlyingtheHDMRrepresentationthatonlylow-ordercorrelationsamongtheinputvariablesarelikelyto
havesignificantimpactsupontheoutputsformosthigh-dimensionalcomplexsystems.HDMRexpansionscanbebroadly
categorizedasRS(RandomSampling)-HDMRandCut-HDMR.Inthefirstcase,thecomponentfunctionsaredeterminedby
integrating over randomly scattered sample points. Whereas the component functions in Cut-HDMR expansion are
evaluatedalonglinesorplanesorvolumes(i.e.,cuts)withrespecttoareferencepointinsamplespace.
Recently,authorshavedevelopedtheHDMR-basedfuzzyfiniteelementmethodforuncertaindynamicalsystem[29],
where HDMR was constructed based on regularized samples and conditional responses at selected sample points.
Followingarethedistinctivefeaturesofthepresentapproach:(i)randomlyscatteredsamplepointsarechosenfrominput
space,ratherthanalonglinesorplanesorvolumes(i.e.,cuts)withrespecttoareferencepointinsamplespace,(ii)sample
pointsaregeneratedusingSobol’ssequence[30,31], and(iii)differentanalyticalbasisfunctions,includingorthonormal
polynomials, cubic B splines, and polynomials for approximating the RS-HDMR component functions can be used. The
outlineofthepaperisasfollows.Section2givesabriefbackgroundoffuzzyvariables.Structuraldynamicanalysiswith
fuzzy parametric uncertain variables is discussed in Section 3. The proposed high-dimensional model representation
approachisdescribedinSection4.ThreenumericalexamplesillustratingtheproposedmethodisgiveninSection5.The
method is then integrated with a commercial Finite Element (FE) software. In the numerical examples, triangular
membershipfunctionshavebeenusedandtheresultshavebeenvalidatedagainstdirectMonteCarlosimulations.Based
onthestudiesinthepaper,conclusionsaredrawninSection6.
2. Abriefoverviewoffuzzyvariables
The mathematical description of fuzzy variables can be given in various levels of abstraction. The easiest way to
describeafuzzyvariableisviaintervalvariables.AscalarintervalvariablexIcanbedefinedas
xI¼½xl,xh(cid:2)¼fx2R9xlrxrxhg ð1Þ
The scalar interval variable can be extended to vector- and tensor-valued variables in a natural way. In contrast to the
probabilisticdescription,theentriesofavectorortensorarealwaysconsideredtobeuncorrelated.Therefore,generalizing
(1),an-dimensionalintervalvectorsimplylieswithinan-dimensionalhypercubeinRn.Theconceptofthefuzzysetwas
proposed by Zadeh [7] to represent imprecise information in a mathematically rigorous manner. Since then numerous
booksandpapershavediscussedthedetailoffuzzysets(see,forexample,thereviewpapers[8–10]).
R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17 7
Afuzzysetcanbedefinedbyextendingthedeterministic(crisp)valueandcombiningintervalalgebra.SupposeXisa
(real) space whose elements are denoted by x. The key idea in fuzzy set theory is the membership function mA~ :
X-MD½0;1(cid:2)whichmapsXintothemembershipspaceM.IfmA~ðxÞ¼1,thenxisdefinitelyamemberofthefuzzysetA~.On
theotherhand,ifmA~ðxÞ¼0,thenxisdefinitelynotamemberofthefuzzysetA~.Foranyintermediatecase0omA~ðxÞo0,the
membershipisuncertain.Therefore,usingthemembershipfunction,afuzzysetcanbedefinedas
A~ ¼fðx,mA~ðxÞÞ9x2Xg ð2Þ
Ifsupx2XmA~ðxÞ¼1,thenthefuzzysetA~ iscallednormal.Inthispaperonlynormalizedfuzzysetsareconsidered.
The notion of support and a-cuts of a fuzzy set play a crucial role in the computational and analytical methods
involvingfuzzysets.ThesupportofafuzzysetA~ isthecrispsetofallx2X,suchthatmA~ðxÞ40,ormathematically
suppðA~Þ¼fx9mA~ðxÞ40,x2Xg ð3Þ
Thea-cutora-levelsetofA~ isthecrispsetAa suchthat
Aa¼fx9mA~ðxÞZa,x2X,0oar1g ð4Þ
In thecase of bounded supports, we also define the 0-cut as the largest closed interval A ¼½infðsuppðA~ÞÞ,supðsuppðA~ÞÞ(cid:2).
0
Often a-cuts are considered as intervals of confidence, since in case of convex fuzzy sets, they are closed intervals
associatedwithagradationofconfidencebetween[0,1].
3. Fuzzyuncertaintypropagationinstructuraldynamics
Theequationofmotionofadampedn-degree-of-freedomlinearstructuraldynamicalsystemwithfuzzyparameters
canbeexpressedas
MðxÞq€ðtÞþCðxÞq_ðtÞþKðxÞqðtÞ¼pðtÞ ð5Þ
wherex2RMisafuzzyvector,pðtÞ2Rnistheforcingvector,qðtÞ2RnistheresponsevectorandM2Rn(cid:3)n,C2Rn(cid:3)nand
K2Rn(cid:3)narethemass,dampingandstiffnessmatrices,respectively.Whenfuzzynessinthesystemparameters,boundary
conditionsandgeometryareconsidered,thesystemmatricesbecomematriceswithfuzzyentries.Uncertaintiesinsystem
(5) are completely characterized by the membership functions of the entries of the vector x. The main interests in
structuraldynamicanalysisare(a)quantificationofuncertaintyintheeigenvaluesandeigenvectorsand(b)quantification
of uncertainty in the response q either in the time domain or in the frequency domain. Few researchers [20–22] have
considered uncertainty quantification in the frequency response function of a linear dynamical systems with fuzzy
variables.Thisproblem,ineffect,hastosolvetheuncertaintypropagationproblem,whichcanbeexpressedinageneral
wayas
y¼fðxÞ2Rm ð6Þ
wherefð(cid:4)Þ:RM-Rm isa smoothnonlinearfunctionoftheinputfuzzyvectorx. Thefunctionfð(cid:4)Þcanbethefrequency
responsefunctionofEq.(5)ortheeigensolutions.
The uncertainty propagation of fuzzy variables through smooth nonlinear functions can be performed using interval
algebrabasedapproachorglobaloptimisationbasedapproach.Intheintervalalgebrabasedapproach,afuzzyvariableis
considered as an interval variable for each a-cut and propagated using classical interval arithmetic [5]. Unless the
functions are simple and the number of variables is small, this approach produces a very large bound for the output
variables. Such overestimates may not be useful for practical design decisions. While for the global optimisation based
approach, for each a-cut two optimisationproblems are solved for each of the output quantities.The first optimisation
problemaimstofindtheminimumvalue,whilethesecondoptimisationproblemaimstofindthemaximumvalueofthe
output quantity. By combining results for all a-cuts, one can obtain the fuzzy description of the output quantities. The
advantage of this approach is that the bounds are ‘tight’ and the wealth of tools are available for mathematical
optimisationcanbeemployedinareasonablystraight-forwardmanner.Thedisadvantageis,ofcourse,thecomputational
cost as two optimisation problems need to be solved for every a-cut and for every output quantity. For these reasons,
efficientnumericalmethodsarenecessarytousethisapproach.Inthispapertheglobaloptimisationbasedapproachis
adoptedforuncertaintypropagation.
For linear structural dynamic problems, typically one is interested in dynamic response in the time or frequency
domain.Thefuzzyeigenvalueproblemcorrespondingtodynamicalsystem(5)canbeexpressedas
KðxÞ/ðxÞ¼oðxÞ2MðxÞ/ðxÞ 8j¼1,...,n ð7Þ
j j j
here oðxÞ and /ðxÞ are jth fuzzy natural frequency and mode shape of system. Assuming the system is proportionally
j j
damped,thefuzzydisplacementresponseinthefrequencydomainwithzero-initialconditionscanbeexpressedinterms
offuzzyeigenvaluesandeigenvectorsas
qðx,oÞ¼ Xn /jðxÞTpðoÞ /ðxÞ ð8Þ
(cid:5)o2þ2iozoðxÞþoðxÞ2 j
j¼1 j j j
8 R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17
hereqðoÞandpðoÞare,respectively,theFouriertransformsofqðtÞandpðtÞ.Inthetimedomain,thedisplacementresponse
ofsystem(5)canbeexpressedas
qðtÞ¼ Xn (Z to1ðxÞ/jðxÞTfðtÞe(cid:5)zjojðxÞðt(cid:5)tÞsinðodjðxÞðt(cid:5)tÞÞdt)/jðxÞ ð9Þ
j¼1 0 dj
wherethedampednaturalfrequenciesaregivenby
qffiffiffiffiffiffiffiffiffiffiffiffi
o ðxÞ¼oðxÞ 1(cid:5)z2 ð10Þ
dj j j
For many practical applications, linear or nonlinear transformations of fuzzy eigensolutions and displacement response
canbeofinterest.
Here we formulate the optimisation problem in the time domain. But the same formulation also applies to the
frequency domain response with minor modifications arising due to the fact that the frequency domain response is
complexinnature.
Suppose I is the set of integers and r2I is the number of a-cuts used for each of the input fuzzy variables
xj, j¼1;2,...,m,inthenumericalanalysis.Therefore,0r½a1,a2,...,ar(cid:2)o1arethevaluesofaforwhichwewishtoobtain
our outputs. Suppose fjðxakÞ, j¼1;2,...,M, are the functions which relate the output yjk with the interval variable xak
associatedwithak forsomekrr suchthat
yjkðtÞ¼fjðxak,tÞ, j¼1;2,...,m, k¼1;2,...,r 8t ð11Þ
The fuzzyuncertaintypropagationproblemcanbeformallydefinedasthesolutionof thefollowingsetof optimisation
problems:
yjkminðtÞ¼minðfjðxak,tÞÞ) 8j¼1;2,...,m, k¼1;2,...,r 8t2½0,T(cid:2) ð12Þ
yjkmaxðtÞ¼maxðfjðxak,tÞÞ
where½0,T(cid:2)denotesthetimeintervalofinterest.Intotal,thereare2mroptimisationproblemsofeverytimeinstantused
inthenumericalprocedure.Thecomputationalcostofeachoftheseoptimisationproblemsdependsonthedimensionality
ofthexa ,thatis,M.Sincem,r,t,andMarefixedbytheproblem,therearetwowaystoreducethecomputationalcost,
k
namely(a)toreducethenumberoffunctionevaluationsbyusingsuperioroptimisationtechniquesor(b)toreducethe
costofeachfunctionevaluation,forexample,by‘replacing’theexpensivefunctionbyasurrogatemodel.DeMuncketal.
[32] used a Kriging based approach. Adhikari et al. [29] proposed a basic HDMR approach. Clearly, if the cost of each
functionevaluationcanbereduced,thentheoverallcomputationalcosttosolvetheuncertaintypropagationproblemwill
bereducednomatterwhatoptimisationalgorithmisused.HereweexplorearandomsamplingbasedHDMRapproachto
efficiently generate the surrogate model. The aim is to ‘build’ the surrogate model with the minimum number of
evaluationsoftheoriginalfunction.Oncethesurrogatemodelisobtained,MonteCarlosimulationisusedtoobtainthe
minimumandmaximumofthefunction.
4. High-dimensionalmodelrepresentation(HDMR)
4.1. Background
ThebasicHigh-DimensionalModelRepresentation(HDMR)forsystemsFuzzyuncertainvariablewasproposedbythe
authors [29]. In this section we briefly review the ordered sampling based HDMR and consequently the new randomly
sampledHDMR(RS-HDMR).HDMRofanarbitraryM-dimensionalfunctionfðxÞ,x2RMcanbederivedbypartitioningthe
identityoperatorI,calledIM intheM-dimensionalcaseandalsointhe1Dcasehereafter,withrespecttotheprojectors
P1,P2,...,PM.Thiscanbeexpressedasfollows[33–40]:
M M M M M
Y Y X Y X X Y
IM¼ ðPmþðI1(cid:5)PmÞÞ¼ Pmþ ðI1(cid:5)PmÞ Psþ ðI1(cid:5)PmÞðI1(cid:5)PsÞ Pp
m¼1 m¼1 m¼1 sam m¼1s¼mþ1 pam,s
|fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
1term ðM1Þterms ðM2Þterms
XM Y YM
þ(cid:6)(cid:6)(cid:6)þ Pm ðI1(cid:5)PsÞþ ðI1(cid:5)PmÞ ð13Þ
m¼1 sam m¼1
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}
ðMÞterms 1term
M(cid:5)1
composedof2Mmutuallyorthogonalterms.TheorthogonalrepresentationofEq.(13)isamanifestationoftheHDMRand
canberewrittenas[41,42]
XM X XM
fðxÞ¼f þ f ðxÞþ f ðx,xÞþ(cid:6)(cid:6)(cid:6)þf ðx,x ,...,x Þ¼ ZðxÞ ð14Þ
0 i i ij i j 123...M i 2 M l
i¼1 1riojrM l¼0
R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17 9
where f is a constant term representing the zeroth-order component function or the mean response of any response
0
function fðxÞ. f is the first-order term expressing the effect of variable x acting alone upon the output fðxÞ, and this
i i
functionisgenerallynonlinear.Thefunctionf ðx,xÞisasecond-ordertermwhichdescribesthecooperativeeffectsofthe
ij i j
variablesx andx upontheresponse.Thehigherordertermsgivethecooperativeeffectsofincreasingnumbersofinput
i j
variablesactingtogethertoinfluencetheoutput.Thelasttermf ðx,x ,...,x Þcontainsanyresidualdependenceofall
123...M i 2 M
the input variables locked together in a cooperative way to influence the output. Once all the relevant component
functions in Eq. (14) are determined and suitably represented, then the component functions constitute the HDMR,
thereby replacing the original computationally expensive method of calculating the response by the computationally
efficientmetamodel.Usually,thehigherordertermsinEq.(14)arenegligible[28]suchthatanHDMRwithonlyafewlow
order correlations among the input variables is adequate to describe the output behavior. This in turn results in rapid
convergenceoftheHDMRexpansion[35,36],whichhasafinitenumberoftermsandisalwaysexact[28,39]intheleast-
square senses. Other popular expansions(e.g., polynomial chaos)have been suggested [1], but they commonly have an
infinite number of terms with some specified functions, such as Hermite polynomials [39]. HDMR expansions can be
broadly categorized as RS (Random Sampling)-HDMR and Cut-HDMR. In the first case, the component functions are
determined by integrating over randomly scattered sample points. Whereas the component functions in Cut-HDMR
expansionareevaluatedalonglinesorplanesorvolumes(i.e.,cuts)withrespecttoareferencepointinsamplespace.
4.2. Cut-HDMR
Whenanorderedsamplingtechniqueisusedinordertoevaluatetheresponse,itisknownasCut-HDMRexpansion.To
generate the HDMR approximation of any function, more precisely the cut-center based HDMR, first a reference point
x¼ðx ,x ,...,x Þ has to be defined in the variable space. In the convergence limit, where all correlated functions in
1 2 M
Eq.(14)areconsidered,thecut-HDMRisinvarianttothechoiceofreferencepointx.However,inpractice,thechoiceof
referencepointxisimportantforthecut-HDMR,especiallyifonlythefirstfewterms,sayuptofirst-andsecond-order,in
Eq.(14)areconsidered[43].Sobol[35]showedthatthereferencepointxatthemiddleoftheinputdomainappearstobe
the optimal choice. The expansion functions are determined by evaluating the input–output responses of the system
relativetothedefinedreferencepointalongtheassociatedlines,surfacesorsub-volumesintheinputvariablespace.This
processreducestothefollowingrelationshipforthecomponentfunctionsinEq.(14):
Z
f ¼ dxfðxÞ
0
Z
f ðxÞ¼ dxifðxÞ(cid:5)f
i i 0
Z
f ðx,xÞ¼ dxijfðxÞ(cid:5)f ðxÞ(cid:5)f ðxÞ(cid:5)f ð15Þ
ij i j i i j j 0
whereRdximeanstointegrateoverallMvariablesexceptx andRdxijmeanstointegrateoverallMvariablesexceptx and
i i
x, etc. These integrals are generally evaluated using numerical integration techniques. Substituting the component
j
functionsdefinedinEq.(15)intoEq.(14),thegeneralexpressionofthecut-HDMRcanbeexpressedas
fðxÞ¼ X fðxi1,...,xib;xi1,...,ibÞ(cid:5)ðM(cid:5)bÞ X fðxi1,...,xib(cid:5)1;xi1,...,ib(cid:5)1Þ
1ri1o(cid:6)(cid:6)(cid:6)oibrM 1ri1o(cid:6)(cid:6)(cid:6)oib(cid:5)1rM
þ2ð!MðM(cid:5)(cid:5)bbþ(cid:5)11ÞÞ!! X fðxi1,...,xib(cid:5)2;xi1,...,ib(cid:5)2Þ(cid:5)(cid:6)(cid:6)(cid:6)8ðb(cid:5)1ðÞM!ðM(cid:5)(cid:5)2Þb!(cid:5)1Þ! X fðxi;xiÞ
1ri1o(cid:6)(cid:6)(cid:6)oib(cid:5)2rM 1rirM
7 ðM(cid:5)1Þ! X fðxÞ ð16Þ
b!ðM(cid:5)b(cid:5)1Þ!
1rirM
where b is the order of the cut-HDMR approximation, 1rbrðM(cid:5)1Þ and the þ or (cid:5) sign of the last term in Eq. (16)
corresponds to b being even or odd, respectively. Considering the weak role of the higher-order correlation effects, the
approximationislikelytoconvergeatalowercut-HDMRorder,say,b5M.TheparticularformofEq.(16)forb¼1,2or3
correspondstofirst-,second-orthird-ordercut-HDMRcanbeexplicitlygivenas
f^ðxÞ¼ X fðx;xiÞ(cid:5)ðM(cid:5)1ÞfðxÞ, b¼1 ð17Þ
i
1rirM
f^ðxÞ¼ X fðx,x;xi,jÞ(cid:5)ðM(cid:5)2Þ X fðx;xiÞþ ðM(cid:5)1Þ! X fðxÞ, b¼2 ð18Þ
i j i 2!ðM(cid:5)3Þ!
1riojrM 1rirM 1rirM
f^ðxÞ¼ X fðx,x,x ;xi,j,kÞ(cid:5)ðM(cid:5)3Þ X fðx,x;xi,jÞþ ðM(cid:5)2Þ! X fðx;xiÞþ ðM(cid:5)1Þ! X fðxÞ, b¼3
i j k i j 2!ðM(cid:5)4Þ! i 3!ðM(cid:5)4Þ!
1riojokrM 1riojrM 1rirM 1rirM
ð19Þ
10 R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17
Thetermfðx;xiÞisafunctionofthesinglex component(i.e.,acutalongx throughthereferencepointinthefunction
i i i
space), while the other variables, x (cid:7)x, jai, are fixed at the reference point. In the same manner, fðx,x;xi,jÞ is the
j j i j
observedresponseforallthevariables,x (cid:7)x ,kai,j,fixedatthecutcenterexceptforx andx.Asimilarinterpretation
k k i j
wouldapplytohigherordercut-HDMRterms.
Thenotionoffirst-,second-order,etc.usedintheHDMRdoesnotimplytheterminologycommonlyusedeitherinthe
Taylor series or in the conventional polynomial based approximation formulae. In HDMR-based approximation, these
terminologies are used to define the constant term, or for example, terms withone variable or two variables only. It is
recognized that the lower order (e.g., first-order or second-order) function expansions in the HDMR, do not generally
translatetolinearorquadraticfunctions[36].Eachofthelower-ordertermsintheHDMRissub-dimensional,buttheyare
not necessarily low degree polynomials. The computational savings afforded by the HDMR are easily estimated. If the
HDMR converges at b order with acceptable accuracy and considering s sample points for each variable, then the total
numberofnumericalanalysesneededtodeterminetheHDMRisPb ½M!=k!ðM(cid:5)kÞ!(cid:2)ðs(cid:5)1Þk.
k¼0
4.3. RS-HDMR
Distinct, but formally equivalent, HDMR expansions, all of which have the same structure as Eq. (14), may be
constructed.Whendataareconsideredasrandomlysampled,RS-HDMR[28,34]canbeobtained.ForRS-HDMR,wefirst
rescale the variables xi by some suitable transformations such that 0rxir1 for all i. The output function fðxÞ is then
definedintheunithypercube[28]
KM¼fðx1,x2,...,xMÞ=0rxir1, i¼1;2,...,Mg ð20Þ
ThecomponentfunctionsofRS-HDMRpossessthefollowingforms[28,34]:
Z
f ¼ fðxÞdx
0
KM
Z
f ðxÞ¼ fðxÞdxi(cid:5)f
i i 0
KM(cid:5)1
Z
f ðx,xÞ¼ fðxÞdxij(cid:5)f ðxÞ(cid:5)f ðxÞ(cid:5)f ð21Þ
ij i j i i j j 0
KM(cid:5)2
where dxi and dxij are just the product dx dx ...dx without dxi and dxij, respectively. Finally, the last term
1 2 M
f ðx,x ,...,x Þ is determined from the difference between f(x) and all the other component functions in Eq. (14).
123...M i 2 M
The RS-HDMR component functions satisfy the following condition: the integral of a component function of RS-HDMR
withrespecttoanyofitsownvariablesiszero,i.e.
Z 1
f ðx ,x ,...,x Þdx ¼0, s2fi ,i ,...ig ð22Þ
0 i1i2...il i1 i2 il s 1 2 l
whichdefinestheorthogonalityrelationbetweentwoRS-HDMRcomponentfunctionsas
Z
f ðx ,x ,...,x Þf ðx ,x ,...,x Þdx¼0 fi ,i ,...igafj ,j ,...jg ð23Þ
KM i1i2...il i1 i2 il j1j2...jl j1 j2 jl 1 2 l 1 2 l
Thecomponentfunctionsf ðxÞ,f ðx,xÞ,...aretypicallyprovidednumerically,atdiscretevaluesoftheinputvariablesx,
i i ij i j i
x,yproduced from sampling the output function fðxÞ for employment on the right-hand side (rhs) of Eq. (21). Thus,
j
numericaldatatablescanbeconstructedforthesecomponentfunctions,andtheapproximatevalueoffðxÞforanarbitrary
pointxcanbedeterminedfromthesetablesbyperformingonlylow-dimensionalinterpolationoverf ðxÞ,f ðx,xÞ,... :To
i i ij i j
construct the numerical data tables for the RS-HDMR component functions [28,34], one needs to evaluate the above
integrals. Evaluation of the high-dimensional integrals in the RS-HDMR expansion may be carried out by Monte Carlo
randomsampling.Forinstance,NsamplesoftheM-dimensionalvectorxðsÞ¼ðxðsÞ,xðsÞ,...,xðsÞÞðs¼1;2,...,NÞarerandomly
1 2 M
generateduniformlyinKM,andthenf isapproximatedbytheaveragevalueoffðxÞatallx(s):
0
f ¼Z fðxÞdx(cid:8) 1XN fðxðsÞÞ ð24Þ
0 KM Ns¼1
WhenN-1,anaccuratevalueoff canbeobtained.Veryoftentheintegralsconvergequitefast,andamodestvalueofN
0
may give a very good result. Moreover, the approximation of an integral by Monte Carlo sampling often does not
depend significantly on the dimension M. This property is extremely beneficial for high-dimension systems. Direct
determination of all RS-HDMR component functions [28,34] at different values of x, x,yby Monte Carlo integration
i j
requires a large number of random samples. For example, to determine f ðxÞ, different sets of samples of fðx,xiÞ at
i i i
ðx,xiÞðsÞ¼ðxðsÞ,xðsÞ,...,xðsÞ ,xðsÞ,xðsÞ ,...,xðsÞÞwithdistinctfixedvaluesofx areneeded,i.e.
i 1 2 i(cid:5)1 i iþ1 M i
f ðxÞ¼Z fðxÞdxi(cid:5)f (cid:8) 1XN fðx,xiÞðsÞ(cid:5)1XN fðxðsÞÞ ð25Þ
i i KM(cid:5)1 0 Ns¼1 i Ns¼1
R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17 11
Ifx takesmdistinctvalues,thenmNsamplesarenecessarytoconstructthef ðxÞnumericaltable.Similarly,toconstruct
i i i
the f ðx,xÞ numerical table, different sets of samples of fðx,x,xijÞ at ðx,x,xiÞðsÞ¼ðxðsÞ,xðsÞ,...,xðsÞ ,xðsÞ,xðsÞ ,...,
ij i j i j i j 1 2 i(cid:5)1 i iþ1
xðsÞ ,xðsÞ,xðsÞ ,...,xðsÞÞwithdistinctfixedvaluesofx,x areneeded,i.e.
j(cid:5)1 j jþ1 M i j
f ðx,xÞ¼Z fðxÞdxij(cid:5)f ðxÞ(cid:5)f ðxÞ(cid:5)f (cid:8) 1XN fðx,x,xijÞðsÞ(cid:5)1XN fððx,xiÞðsÞÞ(cid:5)1XN fððx,xjÞðsÞÞ ð26Þ
ij i j KM(cid:5)2 i i j j 0 Ns¼1 i j Ns¼1 i Ns¼1 j
þ1XN fðxðsÞÞ ð27Þ
N
s¼1
If both x and x take m distinct values, then m2N random samples are necessary to construct the f ðx,xÞ. The required
i j ij i j
numberofsamplesincreasesexponentiallywiththeorderoftherequiredRS-HDMRcomponentfunctions.Thus,thedirect
approach is prohibitively expensive for the construction of high-order RS-HDMR component functions. To reduce the
sampling effort, the RS-HDMR component functions are approximated by expansions in terms of a suitable set of basis
functions,suchasorthonormalpolynomials
k
f ðxÞ(cid:8) Xaij ðxÞ
i i r r i
r¼1
f ðx,xÞ(cid:8) Xl Xl0 bij j ðx,xÞ ð28Þ
ij i j pq pq i j
p¼1q¼1
wherek,l,l0 aretheintegers,ai,bij aretheconstantcoefficientstobedetermined,jðxÞandj ðx,xÞaretheone-and
r pq r i pq i j
two-variablebasisfunctions.Withtheseformulas,Eq.(14)canbeapproximatedas
f^ðxÞ¼f þXM Xk aijðxÞþ X Xl Xl0 bij j ðx,xÞþ(cid:6)(cid:6)(cid:6) ð29Þ
0 r r i pq pq i j
i¼1r¼1 1riojrMp¼1q¼1
Wehaveemployedthesolutiontoevaluatethecoefficientoforthonormalpolynomialsfromthereferencegivenin[44].
4.4. Sampling
The generation of sample points is crucial for construction of the RS-HDMR based surrogate model. The quality of
sample points governs the rate of convergence of Monte Carlo simulation, which is used to construct the component
functionsofRS-HDMR[45].Quasi-randomsamplepointsaregeneratedusingthejointprobabilitydistributionoftheinput
parameter [45]. Each sample point is chosen independently of all other sample points. In this work, quasi-random
sequencesareusedtogeneraterandomsamplepoints.Thesearetotallydeterministicandhavinglowdiscrepancy.This
ensures a moreuniform distributionof sample pointsin theinput domaincompared to pseudo-randomsample points,
thusattainingfasterconvergence.Thebestknownquasi-randompointsareHaltonsequence[46],Sobol’ssequence[30,31]
andFauresequence[47].Sobol’ssequence[30,31]exhibitsbetterconvergencethaneitherFaureorHaltonsequence[48].
4.5. Summaryoftheproposedmethodology
Basedonthediscussionsofar,theproposedmethodologyfortheRS-HDMRapproximationiseasilyimplementedusing
thefollowingsteps:
1. IdentifytherangeofvariabilityoftheparametricFuzzyvariables.
2. Generatesamples(N)forfuzzyinputparametersfromSobol’ssequence[30,31].
3. Estimatetheresponsesatallthesesamplepoints.
4. Determinef ,f ðxÞ,f ðx,xÞusingEqs.(24)and(28),respectively,usingorthonormalpolynomialbasisfunctions.
0 i i ij i j
5. Constructfirst-andsecond-orderRS-HDMRapproximation,respectively,usingthefollowingsteps:
f^ðxÞ¼f þXM Xk aijðxÞ ð30Þ
0 r r i
i¼1r¼1
f^ðxÞ¼f þXM Xk aijðxÞþ X Xl Xl0 bij j ðx,xÞ ð31Þ
0 r r i pq pq i j
i¼1r¼1 1riojrMp¼1q¼1
6. PerformMonteCarlosimulationontheapproximatedresponsefunctionf^ðxÞusinguniformvariables.Estimatethe
minimum–maximumrangeoftheresponsefortheparticulara-cut.
7. Repeattheprocedurefordifferenta-cuts.
12 R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17
5. Numericalexamples
TheimplementationoftheRS-HDMRapproachinfuzzyanalysisisillustratedwiththehelpofthreeexamplesinthis
section.Whencomparingcomputationaleffortsinevaluatingtheresponseintervals,thenumberofactualFEanalysesis
chosenastheprimarycomparisontoolbecausethisindirectlyindicatestheCPUtimeusage.ForfullscaleMCS,thenumber
oforiginalFEanalysesissameasthesamplesize.WhentheresponsesareevaluatedthroughfullscaleMCS,theCPUtime
ishighersinceanumberofrepeatedFEanalysisarerequired.However,intheproposedmethodMCSisconductedwithin
the proposed framework. Although the same sample size as for direct MCS is considered, the number of FE analyses is
muchlower.Hence,thecomputationaleffortexpressedintermsofFEcalculationaloneshouldbecarefullyinterpreted.
SincetheRS-HDMRapproximationleadstoanexplicitrepresentationofthesystemresponses,theMCScanbeconducted
foranysamplesize.SamplerequiredforRS-HDMRconstructionaregeneratedfromSobol’ssequence[30].
5.1. Example1:modalanalysisofasimplifiedaircraftwing
This example involves the estimation of bounds of the natural frequencies of a vibrating aircraft wing with length
L ¼1mandwidthL ¼0.5m,asshowninFig.1.Themassoftheengineisattachedtotheplate,whichrepresentsthewing,
y z
throughaspringatone-thirddistancefromtherootofthewing.Poisson’sratiofortheplateistakenas0.33.Themassof
theengineistakenas0.80timesthemassofthewing.Thestiffnessofthespringiscalculatedfromthesecondnatural
frequency of the wing. The thickness, mass density and elastic modulus of the wing are taken as fuzzy variables. The
boundsofthevariablesaretabulatedinTable1.Theresponsequantities,whicharethenaturalfrequenciesofthewing
structure,arerepresentedbyanRS-HDMRexpansion.Notethatinthisexample,theresponses(naturalfrequencies)are
monotonicfunctionsoftheparameters.Furthermore,thereisnononproportionaldamping,modecrossingorveeringin
theparameterrangeconsidered,whereeigensolutionscanbeextremelysensitivetotheparameters[50].Theproposed
approach with RS-HDMR approximation and full scale MCS using the commercial FE code ADINA [49] are employed to
evaluatethefuzzycharacteristicsofthenaturalfrequenciesofthestructure.Therefore,computationalefficiency,evenfor
thissimplifiedmodel,isa majorpracticalrequirementinsolvingthedynamicsystem.Usingsamplesgeneratedforthe
RS-HDMR,approximationerrors(comparedwithfullscaleMCS)intheresponseboundsarepresentedinFig.2.Theerrors
intheresponseboundsata-cutsof0and0.2arequitehigh((cid:8)326%),however,theseerrorscanbereducedbytaking
moreRS-HDMRsamplepoints(increaseN),whichdemandsmoreFEanalyses.Thus,thenumberofsamplepoints(N)for
the RS-HDMR approximation will vary depending on the user’s accuracy requirement. The full scale MCS results are
obtained from 105 FE analyses. In contrast, only 256 (input sample for fuzzy parameters are generated from Sobol’s
sequence[30])FEanalysesarerequiredbyRS-HDMRforeacha-cut.Inthisproblem,foreacha-cut,thetotalCPUusagefor
fullscaleMCSis10001s,whereasforHDMRapproximationitisonly46s.
Fig.1. Thefiniteelementmodelofasimplifiedaircraftwingwiththeenginemodeledasalumpedmassattachedtothewingviaaspringat1/3span.
ThestructureismodeledintheADINA[49]FEsoftware.The10(cid:3)5FEmodeloftheplateconsistsof502-Dsolidelementsand303nodes.Afixed
boundaryconditionisappliedattheleft.
Table1
BoundsoftheinputfuzzyvariablesofExample1ata-cut¼0.
Variable xmin xmax
Young’smodulus(GPa) 29 109
Thickness(m) 0.1 0.5
Density(kg/m3) 2000 4000
R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17 13
7 2.5
Mode 1 Mode 1
Mode 2 Mode 2
6 2
Mode 3 Mode 3
Mode 4 Mode 4
5 Mode 5 1.5 Mode 5
%) 4 %) 1
or ( or (
Err 3 Err 0.5
2 0
1 −0.5
0 −1
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8
α−cut α−cut
Fig.2. ErrorsintheresponseboundsofthenaturalfrequenciesusingproposedRS-HDMRbasedapproachcomparedtoMCS:(a)errorsinthelower
boundsand(b)errorsintheupperbounds.
Fig.3. Thefiniteelementmodelofacantileverboxbeam.ThestructureismodeledintheANSYS[51]FEsoftware.
Table2
Firsteightfrequenciesofthecantileverboxbeam(Fig.3).
1 2 3 4 5 6 7 8
Frequencies(Hz) 0.803 2.485 4.103 8.776 9.752 10.576 10.711 12.519
5.2. Example2:modalanalysisofacantileverboxbeam
Considerthecantileveredboxbeam[25]showninFig.3.Thisstructureiscomposedof5platesandmodelledwith516
shellelements(SHELL63)and3624degreesoffreedominthecommercialFEcodeANSYS[51].Itsmaterialpropertiesare
YoungsmodulusE¼72(cid:3)104N=m2,Poisson’sration¼0:3,anddensityr¼2700kg=m3.Itsoverallgeometriccharacter-
isticsarelengthL¼3m,widthl¼0.2m,andheighth¼0.1m,withtheindividualplatesebeing2mmthick.Thefirsteight
frequenciesof thestructurearepresentedinTable 2,which isinexcellentagreementwithpreviousresults [25]. These
frequenciesconstitutetheinitialdatafortheparametricstudy.Theparametricstudyhasbeenperformedconsidering20%
variationoftheplatethickness.Fivefuzzyparametersaredefined,andthevariationofeightfrequenciesofthestructures
are calculated. Using samples generated for the RS-HDMR, approximation errors (compared with full scale MCS) in the
lowerandupperboundsforthefrequenciesisshowninFig.4.ThefullscaleMCSresultsareobtainedfrom105FEanalyses.
In contrast, only 256 (input sample for fuzzy parameters are generated from Sobol’s sequence [30]) FE analyses are
required by RS-HDMR. The results obtained using HDMR-based approach are very close to the reference data obtained
withZadeh’sExtensionPrinciple[7,25];themaximumerrorinthecalculationsoftheloweranduppermodalquantitiesis
on the order of 0.2%. Fig. 4 shows an excellent agreement of the results obtained using the present approach with the
availableresults[25].Ingeneral,thepercentageerrorissignificantlylowercomparedtoExample1.Thehighererrorfor
thelowervaluesofa-cutinExample1canbeattributedtorelativelyhigherparametricvariabilityasgiveninTable1.
14 R.Chowdhury,S.Adhikari/MechanicalSystemsandSignalProcessing32(2012)5–17
0.2
Lower Bound
0.18 Upper Bound
0.16
0.14
0.12
%)
or ( 0.1
Err 0.08
0.06
0.04
0.02
0
1 2 3 4 5 6 7 8
Frequency
Fig.4. ComparisonofthelowerandupperbounderrorsofestimatedfrequenciesusingfullscaleMCSwith105FEanalysesandpresentapproachfor
a-cut¼0.
1
0.5
0
−0.5
−1
0 0.5 1 1.5 2 2.5 3 3.5 4
t (sec)
Fig.5. Problemstatement(Example3):(a)baseisolatedstructurewithanequipmentisolationsystemonthesecondfloorand(b)accelerationhistory.
5.3. Example3:time-domainresponseofasixdegreesoffreedomsystemdynamicalsystem
Thisexampleconsidersafour-storeybuildingexcitedbyasingleperiodsinusoidalpulseofgroundmotion,studiedby
GavinandYau[52].Fig.5(a)showsthefour-storeybuildingwithisolationsystemsandFig.5(b)presentstheacceleration
history.Thebuildingcontainsisolatedequipmentrestingonthesecondfloor.Themotionofthegroundfloorisresisted
mainly by base isolation bearings [52] and if its displacement exceeds D ð¼0:50mÞ then an additional stiffness force
c
contributes to the resistance. The Mass, stiffness and damping coefficients, m, k and c, respectively, at each floor are
f f f
assumedtobeidentical.
Therearetwoisolatedmasses,representingisolated,shock-sensitiveequipmentrestingonthesecondfloor.Thelarger
mass m ð¼500kgÞ is connected to the floor by a relatively flexible spring, k ð¼2500N=mÞ, and a damper,
1 1
c ð¼350N=m=sÞ, representing the isolation system. The smaller mass m ð¼100kgÞ is connected to the larger mass by
1 2
arelativelystiffspring,k ð¼105N=mÞ,andadamper,c ð¼200N=m=sÞ,representingtheequipmentitself.Descriptionsof
1 2
thefuzzyvariablesarelistedinTable3.Thelimitstate isdefinedbythecombinationofthreefailuremodesleadingto