12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 1 A Dynamic Bayesian Network Framework for Risk Assessment of Systems Based on Sensor Measurements Iris Tien Assistant Professor, School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA, USA Matteo Pozzi Assistant Professor, Department of Civil and Environmental Engineering, Carnegie Mellon University, Pittsburgh, PA, USA Armen Der Kiureghian Professor, Department of Civil and Environmental Engineering, University of California, Berkeley, Berkeley, CA, USA ABSTRACT: In this paper, a framework based the dynamic Bayesian network (DBN) is proposed to dynamically monitor the response of structures to hazards. The methodology enables the probabilistic analysis of the response of a system to a hazard that is stochastic, e.g., an earthquake ground motion, as it dynamically evolves in time, based on sensor measurements that are uncertain. The developed DBN framework is applied to the estimation of the inter-story drift of a multi-story shear-type building under earthquake hazard based on accelerometer measurements. Using a simulation approach, the ability of the method to accurately assess the system response is shown. In addition, robustness of the method to system uncertainties, including uncertainties in the structural characteristics, is demonstrated. 1. INTRODUCTION The analysis of the risk and resilience of a system subject to hazards requires probabilistic modeling of the system as well as assessment of the system through the evolution of the hazard. As the use of sensors becomes increasingly widespread in the monitoring of these systems, the challenge is to use the information from these sensors to perform this system risk assessment. A methodology based on a dynamic Bayesian network (DBN) framework is proposed. Specifically, we are interested in monitoring the response of structures to hazards that cause dynamic effects, such as earthquake ground motions. The objective is to accurately assess the structural response as it evolves through time under an excitation that is stochastic and based on sensor measurements that are uncertain. The DBN enables probabilistic analysis of the dynamically evolving structural system. This is beneficial for performing system risk assessment in real-time, and in the time after the hazard has occurred. In providing real-time estimates of the response, the proposed framework informs decision making in the management of structures subject to seismic hazard, for example, to trigger emergency shut-down procedures for a structure once a response variable such as inter-story drift exceeds a given safe threshold. In post-hazard system assessment, the framework processes information collected from sensors to estimate not only the state of the structure at the end of the hazard event, but to map the trajectory of the evolution of the system state, informing, for example, varying safety procedures. Finally, as the proposed analysis depends on assumed values of the system parameters that are subject to uncertainty, including the structural characteristics, we show 12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 2 the proposed inference method to be robust relative to uncertainties in these parameters. 2. BACKGROUND AND RELATED WORK Several studies have used Bayesian methods in structural health monitoring (SHM) to perform damage detection in a structure, including Yang et al. (2006), Vanik et al. (2000), and Sohn and Law (1997). In addition, Bayesian techniques have been used in SHM to identify the modal parameters of a structure. For example, a Bayesian framework is proposed in Yuen and Katafygiotis (2002) to explicitly treat uncertainties in sensor measurements and modeling assumptions to obtain distributions of the modal parameters, including the most probable values of the parameters and their uncertainties. Both Katafygiotis and Yuen (2001) and Au et al. (2013) use data from ambient vibrations for modal identification to obtain updated distributions of the modal parameters. The system identification and damage detection problems described previously that are addressed in these studies are long-term monitoring problems, however. Our interest is in performing inference on the state of the structural system as it is subjected to a specific stochastic hazard. As the DBN performs probabilistic analysis of a dynamic system, it is ideal for SHM applications, where there is uncertainty in the sensor measurements and the structure is dynamically evolving in time. The DBN methodology proposed below performs inference on the dynamically evolving response of a structure when it is subjected to a stochastic excitation, e.g., an earthquake, based on the uncertain information collected from sensor measurements. By monitoring the structural response as it evolves, we are able to assess the risk of the system due to the hazard, obtaining, e.g., the expected inter-story drift response, and the probability of exceeding a maximum allowable response. It is noted that the present analysis assumes linear structural behavior as well as Gaussianity of both the earthquake and ambient vibration input excitations to allow the use of Gaussian models and Kalman filter (KF). As such, it is appropriate for serviceability earthquakes and operating-basis seismic events. 3. PROPOSED METHODOLOGY 3.1. System Formulation We model the dynamical system as a cascaded system of two sub-systems: a ground and a structural dynamical sub-system as shown in Figure 1. Figure 1: Dynamical system model, consisting of ground and structural sub-systems. The ground dynamical sub-system takes a modulated white-noise input ๐ค(๐ก), representing the motion at the bedrock, and outputs the acceleration ๐?(๐ก) on the ground surface. The structural dynamical sub-system takes ๐?(๐ก) as excitation and produces the structural response ๐ฎ?(๐ก), the vector of nodal displacements relative to the ground. Our interest lies not only in inferring the instantaneous values of ๐ฎ?(๐ก), but also in its peak values over time. The equation describing the motion on the ground surface relative to the bedrock is given by ๐ข? + 2๐?๐?๐ข? + ๐??๐ข? = โ๐ค (1) where ๐? and ๐? define the angular frequency and damping ratio of the ground filter and ๐ค denotes the modulated white-noise acceleration at the bedrock. Written in first-order form with ๐ณ? = ๐ข? ๐ข? ?, Equation (1) becomes ๐ณ? = 0 1โ๐?? โ2๐?๐? ๐ณ? + 0โ1 ๐ค (2) The total acceleration at the surface of the ground, ๐?, is obtained as ๐? = ๐ข? + ๐ค = โ๐?? โ2๐?๐? ๐ณ? (3) The equations of motion for a linear structure subjected to base motion is ๐๐ฎ? + ๐๐ฎ? + ๐๐ฎ? = โ๐๐ข๐? + f (4) Ground dynamical sub-system Structural dynamical sub-system w ag us12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 3 where ๐, ๐, ย and ย K denote the mass, damping, and stiffness matrices, respectively. ๐ข is the influence vector, where each element of ๐ข equals the displacement in the corresponding degree of freedom from the application of a unit ground displacement. f models a random external force vector that represents ambient vibrations, i.e., random external loads acting on the structure in addition to the earthquake base input, adding uncertainty to the system. In first-order form, using ๐ณ?? = [๐ฎ?? ๐ฎ?? ย ], Equation (4) becomes ๐ณ? = ๐ ๐โ๐??๐ โ๐??๐ ๐ณ? + ๐โ๐ข ๐? + ๐๐?? f (5) Combining Equations (2), (3), and (5), we obtain a representation of the full dynamical system in first-order form. Defining the system state vector as ๐ณ? = [๐ณ?? ๐ณ??], we have ๐ณ = 0 1โ๐?? โ2๐?๐? ๐ ย ย ย ย ย ย ย ย ย ย ย ย ๐๐ ย ย ย ย ย ย ย ย ย ย ย ย ๐๐ ๐๐ข๐?? ๐ข2๐?๐? ๐ ๐โ๐??๐ โ๐??๐ ๐ณ+ ย ย ย ย ย ย ย ย 0โ1๐๐ ๐ค +๐๐๐๐?? f (6) Consistent with previous studies (e.g., Gasparini et al. 1983), we define the matrix ๐? and vector ๐? in a state-space representation of the system in continuous time such that ๐ณ = ๐?๐ณ+ ๐?๐ค +๐?f. Discretizing in time domain in the state-space framework requires the standard transformations ๐? = ๐๐?โ? , ๐? = ๐??๐(๐? โ๐)๐? , and ๐? = ๐??๐(๐? โ ๐ฐ)๐? (Johansson et al. 1999). This leads to ๐ณ??? = ๐?๐ณ? + ๐?๐ค? + ๐?๐? (7) as the full equation of motion for the system in discrete time step ๐. For the sake of simplicity, hereafter we drop the subscripts ๐ and write the discretized equation as ๐ณ??? = ๐๐ณ? + ๐๐ค? +๐๐?. 3.2. Modeling the Excitation To represent the non-stationarity of the ground motion, we model the acceleration at the bedrock as a modulated, band-limited white-noise process. Thus, ๐ค(๐ก) is normally distributed with zero mean and a time-varying variance ๐??(๐ก). Following Rezaeian and Der Kiureghian (2010), ๐??(๐ก) is taken as proportional to a gamma probability density function (PDF), with ๐ก = 0 indicating the beginning of the seismic event. The gamma PDF is a reasonable model for this purpose, since it starts at and tends to zero, is zero for negative values, and the shape is skewed with a longer right tail, which is typical of earthquake motions. The parameters of the gamma PDF are determined in terms of descriptive variables of the seismic event. Specifically, we take the mode of the distribution to coincide with the time of the maximum intensity of the ground motion, ๐กโข โข? , and the middle 90% of the distribution to represent the effective duration of the earthquake motion, ๐ท?? โข , which we define as the time between 5% and 95% Arias intensity values. After derivation, these modeling assumptions lead to the shape and scale parameters of the gamma distribution as ๐ = ? โข โข?? + 1 and ๐ = โ ?? ๐กโข โข? + ?? 4๐กโข โข? ? + ๐ท?? โข ? , respectively. The modulating function is scaled by a factor to achieve the desired intensity of the motion. Figure 2 shows an example simulation of the ground acceleration ๐?(๐ก). Figure 2: Ground acceleration ๐? with time-varying variance ๐??(๐ก) proportional to the gamma PDF. 0 10 20 30 40 50โ1.5โ1โ0.500.511.5time [s]ground acceleration [m/s2 ]12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 4 This example excitation is a realization of a modulated white-noise process with earthquake parameters set at ๐กโข โข? = 20 ย seconds ย (s) and ๐ท?? โข = 25 ย s, and multiplied by a scale factor of 40. This results in a root mean square of the ground acceleration, calculated over the duration of the response 0, 50 ย s , of 0.64 ย ๐/๐ ?. For this example, the discretization time step is set as ฮ๐ก = 0.01 ย s. This effectively cuts frequencies off at 50 Hz. 3.3. System Evolution and Observation Equation Given the state-space representation of the system, the system evolution from time step ๐ to ๐ + 1 is described as in Equation (7). We take ๐? to be normally distributed with zero mean and covariance matrix ๐??๐ , with statistically independent values for different time steps, and values independent of the ground motion and initial condition. Including this additional force increases the uncertainty in our model. We assume the sensors mounted on the structure measure total acceleration, which is given by ๐? = ๐ฎ? + ๐ข๐? (8) From Equation (4) and ๐ณ as defined in Equation (6), we obtain ๐? = โ๐??๐๐ฎ? โ๐??๐๐ฎ? +๐??๐ ย ย ย ย ย ย ย ย ย ย ย ย ย ย ย ย ย = โ๐?? ๐ ๐ ๐ ๐ ๐ณ+๐??๐ (9) Let ๐ define a matrix that selects the degrees of freedom, where accelerometers are placed. The observation equation at each time step ๐ is then given by ๐ฒ? = ๐๐ณ? + ๐? + ๐๐??๐? (10) where ๐ = โ๐๐?? ๐ ๐ ๐ ๐ is the observation matrix and ๐? denotes the vector of measurement errors. We take the measurement errors ๐? to be normally distributed with zero mean and time-independent common variances ๐?? , and assume errors at different times and different locations are independent. The objective is to use these sensor measurements to infer the response of the structure as it evolves with time under seismic excitation. 3.4. Dynamic Bayesian network With a stochastic excitation, a dynamically evolving structural response, and uncertainties in the system, including uncertainties in sensor measurements, the system to be analyzed is well modeled as a DBN. Figure 3 shows the representation of the system as a DBN. Figure 3: Representation of the system as a DBN. The DBN consists of a sequence of BNs, each representing the system at a slice in time, 1,โฆ , ๐ . Terms with zero subscripts represent initial values at time zero, which are also uncertain. The DBN shows the evolution of the system state ๐ณ over time, taking ๐ค and ๐ as input. The measurements ๐ฒ are then taken from ๐ณ with measurement error ๐. The observed variables ๐ฆ? are indicated by the shaded nodes. As stated previously, we assume a linear Gaussian system to allow the use of the Kalman filter (KF) to solve the DBN. That is, each of the random variables is assumed to be conditionally Gaussian given its parent(s) with the mean defined by a linear function of the parent(s). For concision, the details of the KF and Kalman smoother (KS) are not discussed here. The reader is referred to sources including Welch and Bishop (1995) and Murphy (2002) for more information. Using the KF to process the information in the DBN, we obtain estimates of the system state. It is these estimates that we use to probabilistically assess the response of the system under seismic load. y1 z0 w0 z1 โฆ โฆ f0 ฮฝ1 yk+1 wk-1 zk wk zk+1 fk ฮฝk+1 yk ฮฝk wn-1 zn yn ฮฝn fk-1 fn-1 12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 5 4. RESULTS 4.1. Example System We apply the proposed method to a 10-story, shear-type building as shown in Figure 4. Figure 4: 10-story shear-type building model. The building is of nominally uniform mass and stiffness with these parameter values set such that the nominal fundamental period is 1 ย s. We assume the building is modally damped with nominal damping ratios being 5% in each mode. The parameters of the ground filter are set at ๐?/2๐ = 1.5 ย Hz and ๐? = 0.4 . We take the stochastic excitation at the bedrock level as described for the example realization in Figure 2, with earthquake parameters ๐กโข โข? = 20 ย s and ๐ท?? โข = 25 ย s. The estimated mean peak ground acceleration is 1.6 ย m/s? , obtained by multiplying the root mean square of the ground acceleration by a peak factor of 2.5. We use simulation to investigate the accuracy of the proposed method. We first simulate a seismic event as well as the ambient noise. We compute the structural response from this generated combined excitation, and we call this the โactualโ response. We then simulate measurements of this response, including measurement noise. These represent the observations of floor accelerations that we obtain from the accelerometers mounted on the structure. Then, assuming we have only these noisy measurements of floor accelerations and the system and ground motion parameters, we use our dynamic Bayesian network (DBN) formulation of the system to estimate the response of the structure to the stochastic seismic loading. The proposed formulation enables us to perform probabilistic inference on the relative displacements of the structure from our estimate of the system state ๐ณ? . In this study, we are interested in the inter-story drift response in particular. While we are able to analyze the inter-story drift throughout the building, for consistency in the results presented in the following sections when looking at time trajectories of the response, we will be examining one inter-story drift in particular, inter-story drift #5 between floors 4 and 5. 4.2. Inference With one sensor placed on the top floor of the building, Figure 5 shows the time history of the inter-story drift #5, including the โactualโ value and the Kalman filter (KF) and Kalman smoother (KS) estimates. The standard deviation of measurement error is set at ๐? = 0.5 ย m/s?. Figure 5: Time history of inter-story drift #5: actual response, KF, and KS estimates. Because the time histories are close, Figure 6 shows the results for one particular peak that m10 m10 u10 0 5 10 15 20 25 30 35 40 45 50โ0.03โ0.02โ0.0100.010.020.030.04time [s]interโstory drift # 5 [m] actualKF ยตKF ยตยฑmKS ยตKS ยตยฑm12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 6 occurs at ๐ก = 19 ย s Any other segment of the time history can similarly be chosen to analyze the results of the estimation. Figure 6: Time history of inter-story drift #5: actual response, KF, and KS estimates (peak at ๐ก = 19 ย s). In Figure 6, for the KF and KS estimates, the bold lines indicate the mean estimates while the thinner lines indicate the mean estimates +/โ one standard deviation. Comparing the actual response with the estimates, we see that, consistent with the theory, employing the KS improves the accuracy of the response estimate compared to using only the KF. Utilizing the additional information of the measurements over the entire time history results in a KS mean estimate that is closer to the actual inter-story drift values and a decreased variance in the estimate. It is noted that the computation time for inference for each time slice of the DBN is on the order of 10?? ย s , with an average of 7.6ร10?? ย s. Compared to a time step of 10?? ย s, the speed of the DBN inference enables real-time assessment of the system state as it is evolving under earthquake load. 4.3. Robustness to Uncertainty In the previous analyses, we assumed that parameters of the structure were known. In reality, these parameters are subject to uncertainty. In this section, we investigate the robustness of the estimation results to this uncertainty in the structural parameters. Figure 7 shows the effect of uncertainty in the story stiffnesses of the structure on the accuracy of the estimation of the inter-story drift #5 response over the full time history. In this analysis, we use the nominal values of the story stiffnesses to obtain the KS estimate of the response. The nominal values represent our best guess of the structural parameters, e.g., based on design drawings. For the โactualโ structure, the story stiffnesses are considered as random variables and vary from the nominal values. We sample the actual stiffness of each story randomly from the lognormal distribution with the mean equal to the nominal value and over a range of coefficients of variation (c.o.v.โs), assuming statistical independence from story to story. It is these randomly sampled values that we use in our simulation of the actual response and measured values. Figure 7 shows the actual response and the KS estimate for increasing values of the c.o.v. of story stiffnesses k. Because the time histories are close, Figure 8 shows the results for the peaks that occur between ๐ก โ 19, 24 ย s . Any other segment of the time history can similarly be chosen to analyze the results of the estimation. Figure 7: Time history of inter-story drift #5: actual response versus KS estimates with varying c.o.v.โs of story stiffnesses 0-20%. 18.95 19 19.05 19.1 19.15 19.2 19.25 19.3 19.350.0060.0080.010.0120.0140.0160.0180.020.0220.0240.026time [s]interโstory drift # 5 [m] actualKF ยตKF ยตยฑmKS ยตKS ยตยฑm0 5 10 15 20 25 30 35 40 45 50โ0.0200.02c.o.v. of k = 0% actualKS0 5 10 15 20 25 30 35 40 45 50โ0.0200.02c.o.v. of k = 5% actualKS0 5 10 15 20 25 30 35 40 45 50โ0.0200.02c.o.v. of k = 10%interโstory drift # 5 [m] actualKS0 5 10 15 20 25 30 35 40 45 50โ0.0200.02c.o.v. of k = 15% actualKS0 5 10 15 20 25 30 35 40 45 50โ0.0200.02c.o.v. of k = 20%time [s] actualKS12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 7 Figure 8: Time history of inter-story drift #5: actual response versus KS estimates with varying c.o.v.โs of story stiffnesses 0-20% (peaks in interval ๐ก =19, 24 ย s). In Figures 7 and 8, the c.o.v.โs of the actual stiffnesses range from 0-20%. Because each level of c.o.v. produces a different trajectory for the actual response, we compare each estimate with the corresponding actual response for that c.o.v. value, as shown in each subplot. We see that, as expected, as the c.o.v. increases, the KS estimates diverge from the actual responses. The effect is slight, however. Even as c.o.v. increases up to 20%, we see close correspondence between the actual and estimated inter-story drift responses. Figure 9 shows the root mean square errors (RMSEโs) of the mean KS estimates of inter-story drift #5 compared to the actual responses as a function of the c.o.v.โs of story stiffnesses. The MSE is computed over the duration of the response 0, 50 ย s as the average squared error between the estimated and actual inter-story drift values, and the RMSE is its square root. For Figure 9, the RMSEโs were computed across the range of c.o.v.โs with a step size of 0.01. The circles indicate the values at the 5 c.o.v. levels presented previously, i.e., for c.o.v. = 0%, 5%, 10%, 15%, and 20%. Figure 9: RMSEโs of KS-estimated inter-story drift #5 with varying c.o.v.โs of story stiffnesses 0-20%. In Figure 9, we see that, as we expect, the RMSE increases with the increasing uncertainty in the story stiffnesses. The increase, however, is gradual and relatively small in magnitude. Specifically, the RMSE increases from 7.6ร10-4 m to 1.4ร10-3 m over c.o.v.โs ranging from 0-20%, with the maximum RMSE an order of magnitude smaller than the peak inter-story drift response. Note that the non-monotonicity of the RMSE is due to random variation, which will be reduced by taking the average of many simulations. The RMSE shown here is the result from one simulation at each value of c.o.v. The reason the inference is robust relative to the uncertainty in the structural stiffnesses is due to the updating that is performed with the DBN formulation of the problem. The variability in the structural parameters affects the estimate of the system state. However, at each time step, we have information from the sensor measurements, and we use this information to update the estimation. This updating occurs at every time step, in this case ฮ๐ก = 0.01 ย s. Thus, the estimate 19 20 21 22 23 2400.010.02c.o.v. of k = 0% actualKS19 20 21 22 23 2400.010.02c.o.v. of k = 5% actualKS19 20 21 22 23 2400.010.02c.o.v. of k = 10%interโstory drift # 5 [m] actualKS19 20 21 22 23 2400.010.02c.o.v. of k = 15% actualKS19 20 21 22 23 2400.010.02c.o.v. of k = 20%time [s] actualKS0 0.05 0.1 0.15 0.20.60.811.21.41.61.822.22.4 x 10โ3RMSE [m]c.o.v. of k12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12 Vancouver, Canada, July 12-15, 2015 8 is quickly corrected by the measurement information before the estimate has time to evolve incorrectly based on the incorrectly assumed nominal values, making the method robust to uncertainty in the structural parameters of the system. 5. CONCLUSION We have developed a methodology based on the DBN framework to analyze the response of a structure to seismic loading based on uncertain sensor information. We have shown the method to accurately estimate the structural response, e.g., inter-story drift, based on acceleration measurements as the building is subjected to a seismic hazard. The ability of the framework to assess the system state in real time informs decision making in the management of structures subject to seismic hazard. In addition, we have shown the robustness of the method to uncertainties in the structural parameters of the system. Future work will be to relax the assumption of a linear Gaussian system to analyze nonlinear structural behavior. This will extend the use of the DBN framework beyond serviceability earthquakes and operating-basis seismic events to perform inference on structural behavior under more extreme hazard events. 6. REFERENCES Au, S. K., Zhang, F. L., and Ni, Y. C., โBayesian operational modal analysis: Theory, computation, practice,โ Computers and Structures, Vol. 126, pp. 3-14, September 2013. Gasparini, D. A. et al., โRandom Vibration of Cascaded Secondary System,โ Proceedings of the 7th SMIRT Conference, Chicago, Illinois, p. 445-451, 1983. Johansson, R., Verhaegen, M., and Chou, C. T., โStochastic Theory of Continuous-Time State-Space Identification,โ IEEE Transactions on Signal Processing, Vol. 47, No. 1, pp. 41-51, January 1999. Katafygiotis, L. S., and Yuen, K. V., โBayesian spectral density approach for modal updating using ambient data,โ Earthquake Engineering and Structural Dynamics, Vol. 30, No. 8, pp. 1103-1123, August 2001. Murphy, K. P., โDynamic Bayesian networks: representation, inference and learning,โ Doctoral Thesis, 2002. Rezaeian, S., and Der Kiureghian, A., โStochastic Modeling and Simulation of Ground Motions for Performance-Based Earthquake Engineering,โ Pacific Earthquake Engineering Research Center (PEER) Report No. 2010/02, University of California, Berkeley, June 2010. Sohn, H., and Law, K. H., โA Bayesian probabilistic approach for structure damage detection,โ Earthquake Engineering and Structural Dynamics, Vol. 26, No. 12, pp. 1259-1281, December 1997. Vanik, M. W., Beck, J. L., and Au, S. K., โBayesian probabilistic approach to structural health monitoring,โ Journal of Engineering Mechanics, Vol. 126, No. 7, pp. 738-745, July 2000. Welch, G., and Bishop, G., โAn introduction to the Kalman filter,โ Department of Computer Science, University of North Carolina at Chapel Hill, Technical Report 95-041, 1995. Yang, J. N., Lin, S., Huang, H., and Zhou, L., โAn adaptive extended Kalman filter for structural damage identification,โ Structural Control and Health Monitoring, Vol. 13, No. 4, pp. 849-867, July/August 2006. Yuen, K. V., and Katafygiotis, L. S., โBayesian modal updating using complete input and incomplete response noisy measurements,โ Journal of Engineering Mechanics, Vol. 128, No. 3, pp. 340-350, March 2002.