PROBABILISTICENGINEERINGMECHANICSELSEVIERProbabilistie Engineering Mechanics 15 (2000) 309-315www.elsevier.com/locate/probengmechCumulativeformationofresponsesurfaceanditsuseinreliabilityanalysisP.K.Dasa.*,Y.ZhengbamNaalhiturendEngeeingUnivylawGlaow20Sland"DepartmentofNavalArchitecture,DalianUniversityofTechnology,Dalian6023,ChinaAbstractAn improved response surface method is proposed and applied to the reliability analysis of a stiffened plated structure.The responsesurfacefunction isformed in acumulativemannerinordertoproperlyaccountforthe secondordereffectsintheresponsesurfacewithacceptablecomputational effortinvolvedintheevaluationofthestatefunction.First,a linearresponsesurfaceisformedbysearchingforthedesignpointbythefirstorderreliabilitymethod.Thevectorprojectiontechniqueis usedtoallocatethesamplingpointsclosetotheresponsesurface.Then,thelinearresponsesurfaceisimprovedbyaddingsquareterms,andthesecondorderreliabilitymethodisemployedtosearchforthedesign point.All the available sampling points except those generated in the very initial stage are used to obtain a well-conditionedsystemmatrixforregression.Lastly.theobtainedresponsesurfaceischecked fortheselectedsamplingpoints.Iftheresponsesurfacefunction is not satisfactory,it is adaptivelyimprovedby adding crossterms and removing someof the second orderterms.In this wayanappropriateincompletesecondorderresponsesurfacecanbeobtainedandusedforreliabilityanalysis.Examplesaregiventodemonstratethe features of the proposed method.2000 Elsevier Science Ltd.All rightsreservedKeywords:Reliability;Stiffenedplates;Responsesurface;Regression statistics; Structural mechanics1Introductiong(x)=O at the design point x,which is the most likelyfailure point.The second order reliability method (SORM)In reliability analysis,theperformance of a structure orobtains theapproximatesolutiontoEq.(1)byexpandingany other system is governed by a statefunction of a vectorg(x)=O at xwith the second order effects (the curvatures)ofbasicrandomvariates(r.v.).Theprobabilityoffailureisconsideredexpressedastheintegral of thejointprobabilitydensityIn the case that the statefunction is not explicitlyknown,function over the failure domain.Let the random vectorFORM/SORMarenotdirectly applicable.Therefore,MontebedenotedbyX=(XX,...X,),wherenisthenumberCarlo simulation (MCS)seems to bea suitable surrogate.of basic r.v.,and thejoint probabilitydensity function byThe advantages of MCS are obvious, but an innate dis.fx(x);then thefailureprobability is expressedbyadvantageof itistheformidablecomputationaleffortforproblemsinvolvinglowreliabilityoffailureortheproblemsfx(a)dr(1)P:that require a considerable amount of computation in each(x)s(sampling cycle.Toovercome thisdisadvantage,numerouswhereg(x)is the statefunction and is defined as follows:variance reduction techniques have been proposed, e.g.importance sampling[1,2],directional simulation[3],g(x)=0:thelimit stateg(x)>0:the safedomain :conditional expectation [4],etc.For certain problemsFORM/SORMcan beeffectively combinedwithMCS(2)g(x)≤0:thefailure domainDespite the above-mentioned methods, there are manypractical problemsthat cannotbehandledby these methods.Inpracticeg(x)may or maynot be explicit,and theanalyt-Recourse must,therefore,be made to othermethodologies.ical solution of Eq.(l)is inmanycasesdifficulttoobtain.whichmaynot be as accurate as the abovemethods but areTherefore,various methods and algorithms have beenneverthelessfeasibleforawiderspectrumofproblems.Theproposed to approach this integration.responsesurfacemethod (RSM)issuchamethod.The first order reliability method (FORM) obtains theTheRSM stemmedfrom experimental design andwasapproximate solution to Eq.(1)by thelinearexpansion oflater introduced into numerical simulation in reliabilityassessment of complex multivariate systems.The basic+ Comresponding author. Tel.: +44-141-330-4563; fax: +44-141-946-idea of the RSM is to approximate the actual state function.62030266-8920/00/S-see front matter2000 ElsevierScience Ltd.All rights reserved.Pll:S0266-8920(99)00030-2
P.K.Das,Y.Zheng/ProbabilisticEngineeringMechanics15(2000)309-3310complexproblem.Nonetheless,the conceptofpositioningwhich may be implicit and/or very time-consuming tosampling points close to theresponse surface has a bearingevaluate,with the so-called response surface functiononthepropertyoftheRSE.(RSF)thatis easiertodeal with.The RSF commonlyIn the presentpaper,theresponsesurfaceisgenerated inatakes the form of polynomials of the basic r.v.Regressionstepwisefashion.It is aimed at problems where it is imprac-is usually performed to determine the RSF by the least-tical to conducta regressionbased on thecentral compositesquares method (LSM).After the response surface hasdesign due to prohibitive computational effort.Firstly,abeenfitted to a sequenceofsampling points,thereliabilitylinear response surface is used to search for the designanalysis can then be carried out.The crux of the RSM is topointbyFORM.Thevector projectionproposedbyKimachieve a good fit of theRSF to the samplingpoints,par-and Na [8] is employed with modifications. Then, theticularly in the neighbourhood of the design point. Centralobtained linearresponse surfaceis improvedbyaddingsquarecomposite design is frequentlyused in the RSM.Many othertermsandSORMisappliedtosearchforthedesignpoint.Allversions of RSMhavealso beenproposed.the available sampling points except those generated in theResponse surfaces generally take a quadratic form.very initial stage are used to obtain a well-conditioned LSMHigherorderpolynomials generally are not used for concepsystem equation.Finally,the obtained response surfaceistual as well as computational reasons.Occasionally thefurther improved and checked alternately.If the RSF iscross-product terms are ignored (Bucher and Bourgund)not satisfactory,it is augmented to include cross termsand the quadratic is incomplete.In other cases (Wong etand certain second order terms areremoved fromtheRSFal.)they areincluded.FORM/SORMmaybe inaccurateNo new sampling points are needed in the checkingwhen used with a response surface approximation to theprocess.Therefore,this process can be regarded aslimit state because the gradient vector may be highlypost-processing.It may be very efficient and accurateinaccurate.for time-consuming problems by exploiting all of theWong [5] proposed an RSF containing the cross terms.available sample points.When thereare manyr.v.,thenumberofcross terms canbeThis paper presents two advantageous points.Thevery large.On the other hand, useful information might befirst is to modify and improve the linear responselost whensquared terms areessential.Bucherand Bourgundsurface function based on the vector projection method[6] proposedafastandefficientresponsesurfaceapproachbyKimand Na.Squareterms and cross-producttermsin conjunctionwiththeadvanced MCS.Constant,linearandmay be selected in an iterative fashion, so that bettersquared terms areincludedin theRSF.Only4×n+3accuracy can be obtained while the computational effortsampling points are required to determine the RSF and theis maintained fairly low. The method is capable ofdesign point.Thedesign point is improved by interpolation.treating linear and highly non-linear limit surfaceIt is a natural extension to carry on more iterations inaround one design point (most probablefailurepoint).searching for the design point instead of interpolationThe second and more distinctive point is the cumulative(it maybe an extrapolation [7] as well).However,Kimuse of sampling points.In previous studies on theand Na [8] demonstrated by examples that this algo-response surface method, the samples are discardedrithm produces inaccurate results in some cases, evenafter every step of iteration. This is severe waste ofaftermore iterations.On the contrary,Yang et al [9] foundinformation,soformidable computation is indispensableBucher and Bourgund's algorithm to be very efficient andfor large structural problems.In this paper samples areaccurate.repeatedly used so long as they are close to the currentOne reason why a response surface with squared termscentrepointfor samplingAsmaybeperceived, whenbehavespoorly(Ref.[8])isprobablythattoofewsamplingthe centre point is very close to the design point (i.e.points are used in the formation of the response surface.when the algorithm nearly converges),many samples ofWhen an incomplete polynomial is utilised to fit samplingthe previous iterations contribute morethan do some ofpoints that are just enough to uniquely determine the poly-the current iteration.Economical use of samples makes thenomial,theaccuracyofthefittingwill beverymuchdepen-present algorithm remarkably effective for very time-dent on the locations of the sampling points.Anotherreasonconsuming problems,and this is theobjective ofthe papermaybe that the cross terms are crucial.Awiderange of structural reliabilityproblems may beRecently,Kim and Na [8] proposed a vector projectiontreated by single design point analysis.The treatment ofsampling technique,in which the pointsfor generating themulti-mode failure in reliability analysis is of great impor-linearresponse surface arelocated around thedesign pointtance and is our ensuing research topic.and near to the current response surface. Iteration isperformed to obtain the design point. In their work, anindex representing the non-linearity of the actual state func-2.Searchingtion is introduced to account forthe effects ofnon-linearityHowever,this might not be considered generic because theAs mentioned above, the searching process is performedoptimumf.[8] whichmeasures thedistancefrom thedesignby employing FORM based on a linear response surfacepointtoasamplingpoint,is difficulttoobtainforarealistic
311P.K.Das,Y.Zheng/Probabilistic EngineeringMechanics15(2000)309-315where eg is a small number,then terminate the searchingprocess.Otherwise, setunction(12)h-Vhs-i(3)s-s+1,(x)=a+Zbmand go to Step (2).The initial values of coefficients a,br.b2,..bn are deter-For practical purposes,e'g need not bevery small,nined by LSM using sampling points centred at the initialbecause it does not dictate the final accuracy. In the litera-central pointxand located in eachdirection atture h is generally set to a fixed numberbetween 1 and 3.Inthis paper it is reduced with each searching cycle by taking(4)i=l,2.,nthe square root.This is for the sampling points to spreadx±3/ngxover a range about the central point.Different initial valueThe initial central point can be the mean value point.It mayand the powermay as well be used.For example,a biggeralso be shifted towards the failure domain to expeditepower will be appropriate if one expects more iterations inconvergence.Fora very time-consuming calculation,thethe searching process,and vice versa.sampling points may be one-sided; that is, n +1 pointsare sampled to determine the parameters.FORM is performed to obtainβand the corresponding3.Improvementdesign pointx'byusing the initial RSF.Then the searchingThe improvement process starts from forming an RSFprocess is performed as follows:withconstant,linearand squareterms(5)h-V31. 51, 72.For the sth searching, do until Step (5).(13)Hg"()=a+≥b+Zcr2Ag合EagSag'ag"(bi.b2....bh)a=[AglThis function has n ×2 + 1 coefficients to be determined.Vg=.+."axnax"ax2"The initial g"(x)is formed by using all the sampling points(6)in the searching process except those in the initial stage.There are at least n ×2+1 sampling available,so the3.For each random variateX,i=I,2...,n,do until StepLSM can normallyproceed.If the systemmatrixhappensto sufferrank insufficiency,it is perturbed a littleto make it(4).(7)solvable.This rank insufficiency is to be made up in thewith e,=0e'=(ei,e2,eh)ensuing process.After g"(x)is determined,the SORM iswhere Oy is the Kronecker delta.employed to compute the reliability indexβand thecorresponding design point xThen the improvementT'=e-(e-a)aprocess is carried out as follows:(8)(14)if [T|*Othen=T'T'elsef=eeh—0.9661. 1.2.For the sth searching, do until Step (5).with=[1-(1-8）(9)Vg=(bi+2cixi,b2+2c2x2..bnp=(pi.p...ph)Intheabovequations,nosummationisperformed.i+2cnxn, evaluated at xa small number that perturbs the sampling points slightlyfromtheresponsesurface.Notethattoosmallanemay(15)Q=Vg/Vglresult in ll-conditioned system matrix in LSM.4.Apairofsampling points,X+arelocated inthe ithdirec-3.For each random variate X.i =I,2....n, do until Steption and on either side of the design point x.(4).(10)(16)withe=jX=x±hVnorpe' (ei,e..,eh)wherehng,representsthedistancefrom thesamplingT=e-(e.a)apoints to the design point.With x included, there aretotallynx2+1points,which suffice thedetermination(17)if|T|+0thent=T/Telset=eof the n +1 parameters in g'(x)5.Determine g(x)by use of LSM, and then compute thewithp=[-e(1)]p.= (pi.p...h)reliability index β'and the corresponding design point x*(18)byFORM.Ifwhere e', is similar to e, in the searching process.(11)IBr-B-≤e
312P.K.Das,Y.Zheng/ProbabilisticEngineeringMechanics/5(2000)309-315Table1selected sampling points before being checked.If thisReliability index β of Example1response surface function is unsatisfactory,it willbeFunctionExpressionβbyFORMmodified and checked alternately.This process is as follows.βbySORMg'(c)9.841 + 1.655.x1 +1.8261.Takethe lastdesign point in the improvement process as1.8264.882x2-1.572x3the centre to define a"cube":g"(x)9.081 + 2.368x, +1.8711.9983.041xz-0.162xg +0.232--0.747xz+x-hVnox,≤x≤x,+hVnoxi=1,2....n0.901xg(x)7.000+1.000x+(22)1.9192.0661.000x2+0.999x+1.049×10~5X-where h is user-defined.Count the number of all0.999x2+0.999xsampling points (including those points that belong to0.99x1x2+0.999x1.xthe very initial stage) situated inside the cube. Theg(x)Eq, (25)1.9192.066number of the“in-cube"points is represented by m.2.Fitg"(x)with all the in-cube sampling points.In case thatthe system equation is not well conditioned,chooseabigger h and go back to Step (1).4.Apairof samplingpoints are located in the ith direction3. g"(x) is augmented to form g(x),which has cross terms.and on either side of the design point.X=x±hynoxp80=a+b%+0%+R(19)dixiai(23)i创iThe design pointx is also taken as a sampling point.5.Determineg"(x)using all the available sampling pointsThe initial e is chosen to meetexcept those in the very first stage.Then compute β’andthe corresponding design point xby SORM. If thee≤m-nx2-1(24)system equation is rank sufficientandThecross terms are chosen bythe user.Forexample,theB-r-≤evariates with larger coefficients of variation and the(20)variates having higher correlation may be consideredfirst. The final inclusion of these terms is dependent onwhere egis a small number,then terminate the improve-theregression,Inthis case,the LSM system equation ismentprocess.Otherwise,setpossibly rank deficient. Therefore, the LSM using sin-gular decomposition is employed. The second orderh-(h-1)2s-s+1,(21)terms with coefficients smaller than a given value e arethen removed from g(x) and the ftting is repeated untiland go to Step (2).the system equation is not singular.Like in the searching process, h' is reduced with iter-4.g(x) is checked.For allthe in-cube sampling points xations;differentinitial valueand the powermaybeutilisedi = 1,2....m, if g(x,) bears the same sign as the calcu-inaccordancewiththe expected (orthemaximum allow-latedvalueg(x),forallthepointsg(x)is consideredasable)time of iterations.satisfactory.The augmentation of theRSF, the removalof some of the second orderterms with smallest coef-ficients, and the checking of signs are repeated until asatisfyingRSFhas been accomplished.4. Checking5.Compute the reliability index β and the failure probabil-ity by SORM.Perform sensitivity analysis to clarify theIn the checking process the response surface functionimportance of each parameter in g(x).derived in the previous section is first modified usingItshouldbenotedthat,ifthenumbrofsamplingpointsitoosmallcomparedtothenumberofcoefficients ing(x),orTable 2the actual limit state function is not smooth enough,theDerivatives of RSF w.r.t. xchecking of signs may be unsuccessful. In this case, theRSFAt mean pointAt design pointproposed method will yield a less accurate result.Altern-atively,the weighted LSM may be utilised in order tog'lc)(1.655,4.882,1.572)(1.655, 4.882, 1.572)"bend"the RSF towards the limit state function.In thisg"()(2.368, 3.041, 0.162)(2.070, 5.666, 0.103)g(x)(1, 1, 1)case, the selected sampling points are graded according to(2.686,5.293, 0.0796)thedistance level of a point to the design point.Then the
P.K.Das,Y.Zheng /Probabilistic Engineering Mechanics15 (2000)309-315313Table3Reliability indexβand derivatives of Example2FunctionβbyFORMBbySORMDerivative at design pointDerivative at mean pointg(a)3.42733.4273(0.580, 1.01)(0.580, 1.01)8(n)3.35593.3871(0.593,0.991)(0.782,0.977)g(x)3.35043.3835(0.588,1.00)(0.748,0.957)g(x)3.34973.3825(0.580,1.00)(0.811,1.00)error function is weighted in inverse proportion to thethe mean point and the design point. From the design pointdistance levelsof view, sensitivity at the mean point is of importance.Neither g(x)nor g"(x)gives good approximates; stillg"(x) behaves better than g'(x).5.ExamplesIt should be pointed out that, if more iterations areperformed,g'(x)will produce moreaccurate results,butExample1.Assume a limit state functionthe convergence rate is slow.It is not desirable to carry ong(x,x2,x3)=7++x2+-+赔一x2+xx3more iterations with g'(x) when there have been enoughsamplesforform g"(x) andg(x).(25)in which the random variates all have standard normalExample2.Consider a limit state function in standarddistributions.normal variates as in the example given in Ref. [8].g(xj,x2)= exp(0.2x) + 1.4) - x2(26)Five iterations were performed to obtain the linearresponse surface function g'(x).In order to make compar-ison,the samenumberofsamplingpoints wasmaintainedinThree iterations were performed to obtain the lineartheformation of g"(x)that has square terms and g(x)thathasresponse surface function g(x).The same number ofsquareandmixed terms.The results obtained byusing thesesampling points was maintained in the formation of g"(x)RSFs are listed inTable l,where SORM is based onpoint-and g(x).The results obtained by using these RSFs are listedfitting by improved Breitung's method [10].The results byin Table 3. The results by using the LSF g(x) direct is alsousing the LSF g(x) direct is also included in this table.included in this table.Directional simulation withMonte Carlo simulation with1000000samples yieldsβ=1000000 samplesyieldsβ=3.3834.Because thelimit2.038.It can be observed that the linear RSF is less accuratestate function is not strongly non-linear,all methods yieldbut the results are still comparable with those obtained byacceptable results. It is revealed that g"(x)performs fairlyFORM using LSF.The results of g (x)are very close towell in this case.those calculated by using the LSF.The RSF g(x) is almostthesameas LSFg(x),andhenceproduced thesameresultsExample3.Consider a structural problem of Ref. [11]as g(x).Ifthemixed term x2gis included in theRSF,itsThe limit state function iscoefficient isoftheorderof 10-5XxXsXiIn this particular situation,g(x) is basically the same asg(X)=X2X3X4(27)the LSF.Therefore, it may be used to predict the perfor-XX,mance of the system awayfrom the limit state surface.Itmay also be used for sensitivity analysis.On the otherhand,All the random variables are normal and mutuallyg()andg(x)onlyapproximate theLSFatthedesignpoint,independent.The statistics are listed in Table 4.Threeso they arenot appropriate forpredictions.Table2 showsiterations wereperformed.The reliability indices computedthe derivativeofg'(x)and g"(x)with respect to basicr.v.atvia linear response surfacefunction(RSF)g',RSF withTable4square terms g"and RSF with cross-product terms g areStatistics of Example3given in Table5.The reliability indices computed by directVariableMean valuesStandard deviationTable 5Results of Example 340.010.003X20.300.015EFunctionβbyFORMβBySORM360XX3610226×10~611.3×10~60.417×10~43.4113.411X0.500.050.633×10~53.4153.400X0.120.0060.539×10~63.4133.400X4060.0Eq.(27)3.4133.399
314P.K.Das,Y.Zheng/ProbabilisticEngineeringMechanics15(2000)309-315Table77.5m3.5m7.5mResults of Exampleproblem(Nisthenumberof samplepoints informing文4the response surfacefunction.β,is thereliability indexobtained byFORM,B2 is the reliability index obtained by SORM, βy is the reliability index4Pobtained by directional simulation with 100 000 samples generated from2-2the respective response surface function)P4NFunctionβiB22一β32P500132.0022.0122.02122-0131.6001.629P1.615131.6241.6371.632220000261.6151.6271.627P261.6561.5501.543268-1.654221.4051.307P22zP5moments of inertia are expressed as li=α;A?(α=α2=E3α3=0.0833,α4=0.2667αs=0.2000)TheYoung'sPmodulus is treated as deterministic.E=2.0× 107kN/m2.enElementtypes are indicated inFig.1.Accordingto therules,Pthe limit state function is defined asEPg(A1,A2,A3,A4,A5,P) =0.096 umax(A1,A2,A3,A4,As,P)3m-(29)P5一m3where umax is the maximum horizontal displacement.-PrThree iterations were performed and 13 finite elementanalyses were performed in each iteration. The first 13samples with the mean value point as the centre pointFig. 1. Plane Frame structure.were discarded after the first iteration. The results aregiven in Table 7.Ref. [11] gives β=1.439 by importancesampling.use of Eq.(27)are alsogiven in thetable.Directional simu-lation with100000 samples givesβ=3.4000.Observe thatthe reliability indices differ very little, so it maybe inferred6. Concluding remarksthat in the vicinity of design point the limit statefunction ismoderately non-linear. The root-of-mean-square error EA method for the cumulativeformation of theresponsedecreasesbyan order from g',g"tog.For example,Eforsurface function is proposed. This method includes threetheRSFg',isdefinedassteps,ie.searching,improvementand checking.Itis intendedto tackle the problems that are very time consuming./Theproposedmethod canform aresponse surfacefunc-E=Z(g(x) - g(x)2(28)tion that has constant, linear,and selected quadratic terms.NThe selection of second order terms is based on the leastsquare regression using singular decomposition. ThisExample 4.A3-bay and 12-story frame is idealised withmethod is particularly suitable for theformation of responsebeam elements.Different cross sectional areas A,andsurfacefunction wherethe cross terms are significant.horizontal load P are treated as independent random vari-Insearching and improvementsprocesses,a parameterhables; their statistics are listed in Table 6.The sectionalwhich determines the distance from the central point to asampling point, is reduced with iterations. In this way theTable 6Statistics of Example problemsampling points can spread over a range close to the linearresponse surface so that they later can be used in theVariableUnitMean valueStandard deviationDistributionimprovementand checkingprocesses.Ar0.25Thechecking process is performed based on the sampling0.025lognormalA20.160.016points already obtained and can be viewed as post-proces-lognormalAl0.360.036lognormalsor.The second orderterms are selected adaptivelyto formm2At0.200.020lognormalaproperincomplete second order RSF.Then,forall theAsm20.150.015lognormalselected samplingpoints the signs of theRSF and the actualPKN30.07.5Type I largestperformance function are checked.If the satisfactory RSF
315P.K.Das, Y.Zheng / Probabilistic Engineering Mechanics 15 (2000)309-315[4] Ayyub BM, Chia CY.Generalised conditional expectation for struc-cannot be obtained,the weighted least square method maytural reliability assessment. Structural Safety 1992;11:131-46.be used. Furthermore, the variance analysis can also be[5] Wong FA.Slope reliability and response surface method.Journal ofinvolved in this process.Geotechnical Engincering,ASCE 1985;111(1):32-53,It is also possibleto include inversepowerterms in the[6] Bucher CG, Bourgund U. A fast and eficient response surfaceprocess ofaugmentation.If this is thecase,the SORMmayapproach for structural reliability problems.StructuralSafety1990;7:5766.fail to fare well; therefore,the Monte Carlo Simulation with[7]RajashekharMR,Ellingwood BR.Anew look at the response surfaceVariance reduction techniques should be used.approach for structural reliability analysis. Structural Safety1993;12:20520.[8] Kim S-H, Na SW. Response surface method using vector projectedReferencessampling points. Structural Safety 1997;1:3-19.[9] Yang YS,Lee JO, Kim BJ,Structural reliability analysis usingcommercial FEM package,Proceedings of6th International Offshore[] Bucher CG. Adaptive sampling—an iterative fast Monte Carlo proce-and Polar Engineering Conference,USA,1996.p.3874.dure, Structural Safety 1988:5:119-26.[10] Liu P-L, Lin H-Z, Der Kiureghian A. CALREL User Manual. Report[2] Zheng Y, Fujimoto Y, Iwata M. An empirical fiting-adaptiveNo.UCB/SEMM-89/18.Departmentof CivilEngineering,Universityapproach to importance sampling in reliability analysis. Journal ofof Califomia at Berkeley,CA,USA,1989.Society of Naval Architects of Japan 1991;171:433-41.[11] Zhao G. Reliability theory and its applications for engineering struc-[3] Bjerager P. Probability integration by directional simulation. Journaltures. Dalian: Dalian University of Technology Press, 1996.EngineeringMechanics,ASCE1988:114(8):1285-302