Rare event simulation in complex engineering syste

1RareEventSimulationinComplexEngineeringSystemsSiu-KuiAuInstituteforRisk&UncertaintyandCenterforEngineeringDynamicsSchoolofEngineeringUniversityofLiverpool,UK2RareeventsWhat?Why?How?„Whatisrareevent?„‘Small’‘probability’ofoccurrence„Whybother?„Highconsequence,e.g.,Expectedloss=losswhenfailxfailureprobability„How?„Ignore„Eliminate„Manage3Uncertainty?„Eliminate?=identificationproblem„Model?=riskassessmentproblem4DeterministicAnalysisE.g.,•loading•materialproperties•dimensions•dampingE.g.,•Strain•Stress•Failuretime•TemperatureE.g.,•Finiteelementmodels•Analyticalformulas•EmpiricalformulasResponse)(XYParametersXSystemAnalysis5PropagationofuncertaintyBasedonaprobabilisticapproachSystemAnalysis)(PDFxp)(1XY)(2XY)(3XYMM1X2X3X)(ypYPDF??Response)(XYParametersX6ReliabilityanalysisHowoftendoesthesystem‘fails’?∫==xdxpxIbXYPFPF)()())(()()(PDFypYySpecifyb??„Failureprobability)(xpxbxYF)(:⎩⎨⎧=otherwise0)(if1)(bxYxIFIndicatorfunction7ComputationalDifficulties„Problemcomplexity„Alargeno.ofrandomvariablesX„RelationshipbetweenresponseYandvariableXisimplicit&complicated„Rarefailureevent„Failureshouldbeexceptionratherthanrule8Differentperspectives„Asanintegral„Numericalintegration=notfeasibleinhighdimensions„Asanintegralbutfocusonmaincontributions„Lookfor‘importantregions’andapproximatebasedoninfothere„Gaussianapproximation(similartoLaplacetypeintegrals)„Firstorder/Secondorderreliabilitymethods(FORM/SORM)„Asanexpectation„MonteCarlosimulations„Estimatebyaveragingoverrandomsamples9DirectMonteCarlo10DirectMonteCarlo∫==xdxpxIbYPPFF)()()(∑=≈NkkFXIN1)(1)(toaccordingddistributei.i.d.,...,,21xpXXXN⎩⎨⎧=otherwise0)(if1)(bxYxIFIndicatorfunctionsamplesofno.TotalsamplesfailedofNo.=)(xpxbxYF)(:11DirectMonteCarloPropertyofestimator„Unbiased„Resultsofdifferentrunsscatteraroundtheexactanswer„Variance„Almostsureconvergence„Anyparticularsimulationhistorywillconvergetotheexactanswer∑==NkkFFXINP1)(1~FFPPE=]~[1)~lim(==∞→FFNPPPNPPPFFF)1(]~var[−=12010203040506070809010000.20.40.60.81NJ~(a)DirectMonteCarloExample(PF=0.5)010203040506070809010000.20.40.60.81No.ofsamplesNNJ~(b)FP~FP~Aparticularsim.history100sim.histories13CCDFPerspective[CCDF=ComplementaryCumulativeDistributionFunction]yP(Yy)Np111−=]1[YNp212−=]2[YMNpN11=−]1[−NY0=Np][NYLO1„AsinglerunallowsonetoestimatethewholeCCDFofY„Samples{Y1,Y2,...,YN}=Orderedsamples{Y[1],Y[2],...,Y[N]}14CCDFvsSimulationhistory-4-3-2-10123410-310-210-1100bP(Yb)0200400600800100000.050.10.150.20.25b=1FP~0200400600800100000.020.040.060.080.1b=2FP~0200400600800100001234x10-3No.ofsamplesNb=3FP~15DirectMonteCarloRareevents„Goodfornot-so-rareevents,notefficientforrareevents„Onaveragerequires10failedsamplestohavearelativeerrorof30%infailureprobabilityestimate„E.g.,PF=0.001=onaverage1/0.001=1000samplestoget1failedsample=10x1000=10,000samplesintotal„Canweavoidsimulationofrareevents?NPNPPPPFFFFF1~)1(]~var[2−==δc.o.v.ofestimateFPsmallfor16Importancesampling17Response))((forEstimateyXYP≥10ybDirectMonteCarloSimulationUncertainParameterSpaceXSystemanalysisbXY≥)(18Response))((forEstimateyXYP≥10ybSystemanalysisImportancesamplingsimulationUncertainParameterSpaceX19ImportancesamplingsimulationTheory∫=xdxpxIPFF)()(∫×=xdxfxfxpxIF)()()()(∑=≈NkkkkFXfXpXIN1)()()(1)(toaccordingddistributei.i.d.,...,,21xfXXXNspecified)-(userdensitysamplingimportance:)(xfForproperlychosenf(x),theestimatorhasthesamepropertiesasDirectMonteCarloestimator20Example(ISveryefficient)050010001500200000.511.52x10-3No.ofsamplesNP(Yb)ImportancesamplingDirectMonteCarlo„Importancesamplingisveryefficientwhenitworks!21ImportancesamplingsimulationComments„Efficiencydependscriticallyonthechoiceofimportancesamplingdensity(ISD)„AgoodchoiceofISDrequires„goodunderstandingofproblem=nottrivialforcomplexproblems„Tradeofffailurerateofsamplesandvariabilityofp/f„ImproperchoiceofISDcanleadtobias„ImproperchoiceofISDcansufferfromcurseofdimension)(xf22Bias-example2designpoints„EstimateP(Y3)1x2x}{11bYF=}{22bYF=2323−0*1x*2x33Y=max{Y1,Y2}Y1=(X1+X2)/√2Y2=(-X1+X2)/√2X1,X2i.i.d.N(0,1)23Bias-exampleISDcenteredatdesignpoint1=bias(a)05001000150020000.511.522.53x10-3No.ofsamplesNP(Yb)(b)F1F2∑=≈NkkkkFFXfXpXINP1)()()(1)()(*1xXXfkk−=φ(StandardGaussiancenteredatx1*)24Bias-exampleISDcenteredattwodesignpoints=Nobias0500100015002000012345x10-3No.ofsamplesNP(Yb)(a)(b)F1F2∑=≈NkkkkFFXfXpXINP1)()()(1)(21)(21)(*2*1xXxXXfkkk−+−=φφ(MixtureofstandardGaussianscenteredatx1*andx2*)25}3{=YF0*x3CurseofdimensionExample„EstimateP(Y3)Y=(X1+...+Xn)/√nX1,...,Xni.i.d.N(0,1)∑=≈NkkkkFFXfXpXINP1)()()(1⎥⎦⎤⎢⎣⎡−=∑=−niiinxXssXf12*22/)(21exp)2()(π(Gaussiancenteredatx*withvariances2inalldirections)p(x)f(x)26CurseofdimensionExample010002000012x10-3s=1n=10010002000012x10-3s=0.8010002000012x10-3s=1.2010002000012x10-3n=100010002000012x10-3010002000012x10-3010002000012x10-3n=1000No.ofsamplesN010002000012x10-3No.ofsamplesN010002000012x10-3No.ofsamplesNDimensionnSpreadofISD,ss=1s=0.8s=1.2n=10n=100n=1000∑=≈NkkkkFFXfXpXINP1)()()(1OKOKOKOKSomebiasBreaksdownBreaksdownSomebiasOK27CurseofdimensionComments„Asnincreases,procedureb