Quadratic partial eigenvalue assignment in large-scale stochastic dynamic systems for resilient and economic design

Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Quadratic partial eigenvalue assignment in large-scale stochastic dynamic systems for resilient and economic design S. Das a,n, K. Goswami a, B.N. Datta b a b

Department of Mechanical and Aerospace Engineering, University at Buffalo, Buffalo, NY 14260, USA Department of Mathematical Sciences, Northern Illinois University, DeKalb, IL 60115, USA

a r t i c l e i n f o

abstract

Article history: Received 15 June 2015 Received in revised form 9 September 2015 Accepted 2 October 2015 Available online 31 October 2015

Failure of structural systems under dynamic loading can be prevented via active vibration control which shifts the damped natural frequencies of the systems away from the dominant range of a loading spectrum. The damped natural frequencies and the dynamic load typically show significant variations in practice. A computationally efficient methodology based on quadratic partial eigenvalue assignment technique and optimization under uncertainty has been formulated in the present work that will rigorously account for these variations and result in economic and resilient design of structures. A novel scheme based on hierarchical clustering and importance sampling is also developed in this work for accurate and efficient estimation of probability of failure to guarantee the desired resilience level of the designed system. Finally the most robust set of feedback matrices is selected from the set of probabilistically characterized optimal closed-loop system to implement the new methodology for design of active controlled structures. Numerical examples are presented to illustrate the proposed methodology. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Quadratic partial eigenvalue assignment Structural vibrations Optimization under uncertainty Probability of failure Disjointed failure domains Importance sampling

1. Introduction Infrastructures (e.g., high-rise buildings, bridges, wind-turbines), automotive and mechanical devices (e.g., aircrafts, hard disk drives), etc., are subjected to various types of dynamic loads during their service life. They undergo vibration with high amplitude (resonance phenomenon) when some of their damped natural frequencies are close to the dominant frequencies of the dynamic loads. If the design process for such systems fails to properly account for the effects of the dynamic loads on the response of the systems, the designed structure may collapse (e.g., failure of the Tacoma Narrows bridge, Washington in 1940 [1,2]) or lose functionality (e.g., wobbling of the Millennium bridge, London in 2000 [3,4]). The catastrophic consequences can be prevented through deployment of appropriate vibration control strategies, among which the most popular are the active vibration control (AVC) and the passive vibration control (PVC). The AVC is, however, advantageous over the PVC due to its ability to reduce the vibration level of a structure by modifying only those damped natural frequencies of the structure that lie within the dominant spectrum of the loads. A schematic of the AVC technique and various components of an AVC system employed on a high-rise building are depicted in Fig. 1.

n

Corresponding author. E-mail addresses: [email protected] (S. Das), [email protected] (K. Goswami), [email protected] (B.N. Datta).

http://dx.doi.org/10.1016/j.ymssp.2015.10.001 0888-3270/& 2015 Elsevier Ltd. All rights reserved.

360

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Fig. 1. Schematic of (a) AVC technique and (b) various components of an AVC system.

To understand the working principle of AVC technique, let us consider the governing equation representing the dynamics of structures given by € þCqðtÞ _ þ KqðtÞ ¼ fðtÞ MqðtÞ

ð1Þ

Eq. (1) is obtained via the finite element (FE) discretization of the governing partial differential equations (PDEs) approximating the dynamic response of the physical structural system. In Eq. (1), M; C; K denote the mass, damping, and stiffness matrices of the system, respectively. Also, qðtÞ denotes the displacement response of the system, fðtÞ denotes the dynamic load acting on the system and the dot ðÞ_ denotes time derivative. Note that M; C; K A Mnþ ðRÞ, the set of all real symmetric positive-definite matrices of size n  n where n is the degrees of freedom (dof) associated with the FE model of the system. Since M; C; K are generated by FE technique, they are very large and often inherit nice structural properties, such as definiteness, bandness, and sparsity, which are useful for large-scale computing. The system is, in general, denoted by ðM; C; KÞ and is known as open-loop system. The dynamics of ðM; C; KÞ depend on their damped natural frequencies ωdn and modeshapes xdn . The pair ðωdn ; xdn Þ of the system can be found by solving quadratic eigenvalue problem (QEP) associated with the
 2 quadratic pencil Po ðλÞ ¼ Mλ þCλ þ K , where λ denotes the eigenvalue of the QEP given by Po ðλÞx ¼ 0. Here x denotes the right eigenvector of the QEP and is associated with the eigenvalue λ. Note that λ's occur in complex conjugate pairs and ωdn ¼ ImðλÞ. The mathematical basis of AVC is an inverse eigenvalue problem for the pencil Po ðλÞ, known as the quadratic partial eigenvalue assignment (QPEVA) problem. Suppose that a control force of the form BuðtÞ is applied to the structure for vibration control, where B A Mnm ðRÞ with m r n is a fixed control matrix and the control input vector uðtÞ assumes the form _ þ GT qðtÞ. The matrices F; G A Mnm ðRÞ are feedback gain matrices. Under the application of control force, Eq. (1) uðtÞ ¼ FT qðtÞ then gets modified to

 ð2Þ Mq€ ðt Þ þ C  BFT q_ ðt Þ þ K  BGT qðt Þ ¼ fðtÞ


 The modified system is called a closed-loop system and is denoted by M; C  BFT ; K  BGT . The dynamics of this system are then governed


 2 Mλ þ C  BFT λ þ K  BGT .

by

the

eigenvalues

and

eigenvectors

of

the

closed-loop

pencil

Pc ðλÞ ¼

2p Let the set of eigenvalues of the ðM; C; KÞ be denoted as λi 2n i ¼ 1 with λi i ¼ 1 (with 2p⪡2n) representing the set of problematic 2p ðiÞ eigenvalues which needs to be replaced by μi i ¼ 1 . Let xi ¼ 1 2p be the right eigenvectors of the ðM; C; KÞ associated with the 2p 2n eigenvalues λi 2p i ¼ 1 . Then QPEVA problem is to find feedback matrices F and G such that the spectrum of Pc ðλÞ is fμi i ¼ 1 ; λi i ¼ 2p þ 1 g Þ 2n 2n corresponding to λ remain unchanged. This desired feature is known as no-spillover and the eigenvectors xði i¼ i i ¼ 2p þ 1 2p þ 1 property. In this context, note that, the QPEVA technique is computationally different from the traditional pole/eigenvalue assignment technique for which there exist excellent numerically viable methods [5, Sections 11.2–11.3]. In the past, a handful number of algorithms based on the state-space approach have been developed to implement QPEVA technique for AVC [6–11]. Also, in the last two decades, considerable research efforts have been channelized towards implementation of state-space approach based AVC techniques for buildings [12–17], wind turbines [18–23], hard disk drives [24–26], etc. These state-space approach based AVC techniques, however, turn out to be numerically inefficient for large-scale systems as the matrices associated with the linear model are twice the dimension of M; C; K and they lose the inherent computationally exploitable properties stated earlier. In this regard, the mathematical and engineering challenges are to solve the QPEVA problem in the quadratic setting itself and using only a small number of eigenvalues and eigenvectors that are all that can be computed using the state-of-the-art computational techniques (e.g., Jacobi–Davidson methods [27,28]) or can be measured in a vibration laboratory [29–31]. Furthermore, the no-spillover property must be established by means of mathematical results, because it is not possible to do so computationally or by measurements. Meeting these challenges, several computationally effective algorithms have been developed in the recent years for largescale systems [32–34]. Furthermore, work has been done to compute the feedback matrices such that they have minimum norm (MNQPEVA technique) and the robustness of the closed-loop system under small perturbations of the data is ensured (RQPEVA technique). These algorithms are very well suited in a deterministic setting.

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

361

Uncertainties are, however, invariably induced in the open-loop system matrices M, C, and K from several primitive sources including, but not limited to, variability in input parameters (e.g., material and geometric properties, choice of loading) often modeled as parametric uncertainty [35] and system-level structural complexities typically modeled as model form uncertainty or model uncertainty [36]. In the ensuing discussion, bold notations will be used to indicate random variates (random variables or random vector or random matrix) or deterministic vector/matrix variables. It will be clear from the context whether a bold notation refers to a random variate or a deterministic vector/matrix variable. Deterministic scalar variables will always be denoted by non-bold notations. Realizations or samples of a random variate will be denoted by the corresponding non-bold notations. Variabilities in input parameters may arise due to fluctuation in experimental measurements, scale effect, choice of loading pattern, etc. Structural complexities, on the other hand, are due to variations in structural configurations, topology, imperfect boundary conditions, defective interfaces/joints, etc., at the system or subsystem level. Uncertainties in the openloop system matrices will affect the design of feedback matrices for active controlled systems. The design of feedback matrices will also be affected due to uncertainties in input parameters of control systems. It must be noted that failure to account for the effects of these uncertainties might have adverse effects leading to unstable closed-loop system that will essentially defeat the very purpose of designing active controlled systems to reduce the vibration level. In the recent past, stochastic control theory has been developed to incorporate several sources of uncertainties associated with measured response, external excitation, and the open-loop system into the design of feedback matrices [37]. Few techniques, for example H1 and H2 control, have been devised in this context for design of feedback matrices [37, Sections 5.4 and 7.5]. Over the last two decades, the H 1 =H 2 control techniques have been employed for vibration control of buildings [38–41,17], wind turbines [42], hard disk drives [25,26], etc. A large number algorithms to implement H 1 =H 2 control techniques exist [43–49]. The H 1 =H 2 control algorithms, however, suffer from few drawbacks: (i) they shift all the eigenvalues of the controlled system and (ii) they are computationally inefficient for large-scale systems due to the state space approach adopted. These drawbacks, as compared to the QPEVA technique, are briefly presented in Section 2 using a smallscale and a large-scale example problem. This necessitates the need for development of a computationally effective stochastic QPEVA technique for economic and resilient design of closed-loop systems. The first step towards formulation of a stochastic QPEVA technique is the characterization of uncertainties which can be done efficiently by decomposing the complete system into several subsystems and characterizing the uncertainties associated with each subsystem separately by using the most suitable probabilistic model. For example, conventional parametric probabilistic models [50,35,51] can be used to characterize parametric uncertainties, whereas, recently proposed random matrix theory (RMT) based probabilistic models would be suitable to characterize model uncertainties in subsystems having structural complexities [36,52–54]. Recently, a hybrid probabilistic formalism has been developed to couple the conventional parametric model and the RMT based model [55]. In the present work, we will only focus on parametric uncertainties. In the context of designing resilient and economic closed-loop structural systems plagued by uncertainties, two factors are significantly important: (i) quantification of the desired resilience level in terms of the low probability of failure denoted by Pf and (ii) the cost of implementation of the feedback controllers. This can be achieved by formulating an optimization under uncertainty (OUU) problem [56] that minimizes certain cost function while simultaneously satisfying user-prescribed Pf that depends on application and some desired performance level. Accuracy and efficiency in estimation of Pf is compromised when the failure domain Sf is disjointed and complicated in shape, and not known explicitly as a function of system parameters (for example, available only through large-scale finite element simulations). These difficulties are impossible to address not only using the Monte Carlo simulation (MCS) technique [57,58] but also using the traditional importance sampling (IS) scheme often implemented for computation of low Pf [57–59]. All the above mentioned challenges associated with implementation of AVC within stochastic QPEVA setting are addressed in the present paper through development of a novel technique which consists of three major tasks/contributions: (i) Formulation of methodology for stochastic closed-loop systems: Formulate an OUU problem that will account for various uncertainties associated with a large-scale closed-loop system and result in an economic and resilient system. (ii) Computation of P f : Devise a novel scheme that will identify the unknown, disjointed, and complicated shaped Sf using hierarchical clustering technique [60] and estimate Pf accurately and efficiently based on IS scheme. (iii) Recommendation of design guidelines: Recommend guidelines for choosing a unique set of feedback matrices to ensure robustness in performance of the designed stochastic system. Fig. 2 pictorially summarizes all the tasks and their interdependence. It is worth mentioning that the proposed methodology can also account for the effects of spatio-temporal random loading. In this paper, we will consider only deterministic time-dependent loading. The current paper begins with a comparison of deterministic MNQPEVA and H 1 techniques in Section 2 to illustrate the drawbacks associated with the H 1 technique mentioned earlier and enforce the need for development of a computationally effective stochastic QPEVA technique. Formulation of methodology for economic and resilient design of stochastic closedloop systems based on OUU has been developed in detail in Section 3. The computational difficulty associated with estimation of low probability of failure, which is a part of the OUU problem formulated, has also been identified in Section 3. A novel scheme based on hierarchical clustering technique and importance sampling scheme has been developed in Section 4 that can provide a computationally efficient yet accurate estimate of low probability of failure. The algorithmic details associated with this novel scheme are also presented in this section. Guidelines for choosing a unique set of feedback

362

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Fig. 2. The big picture: summary of all tasks and their interdependence.

matrices for implementation of the stochastic closed-loop system designed are elucidated in Section 5. The methodology formulated in Section 3 is then illustrated in Section 6 through a numerical example with a 4-story steel building and a recorded ground motion. Finally, conclusions are presented in Section 7.

2. Comparison of deterministic MNQPEVA and H 1 techniques In this section, the deterministic MNQPEVA control technique [32–34] is briefly described as it serves as the foundation for the formulation of stochastic QPEVA (see Section 3). A comparison of the deterministic MNQPEVA and H 1 control techniques is also presented here using two example problems, one with a small-scale system and the other with a largescale system, to illustrate the advantages of the MNQPEVA technique over the H 1 technique. 2.1. Brief overview of deterministic MNQPEVA technique A matrix A of size r  s may also be denoted by Ars in this paper for the sake of simplicity. Also, the (i,j)-th element of a matrix A will be denoted by Ai;j in this paper. To present the deterministic MNQPEVA problem, let us assume that the 2p problematic eigenvalue matrix be Λo2p2p ¼ diag λi 2p i ¼ 1 , the target eigenvalue matrix be Λc2p2p ¼ diagfμi gi ¼ 1 , and the matrix 2p of eigenvectors associated with the eigenvalues in λi i ¼ 1 be Xon2p . Using these notations, the parametric expressions for the feedback matrices F and G, as developed in the literature [32], are given by   T ΓTm2p ð3aÞ Fnm Γ ¼ Mnn Xon2p Λo2p2p Z2p2p   T ΓTm2p Gnm Γ ¼ Knn Xon2p Z2p2p

ð3bÞ

where Z is the unique solution of the Sylvester equation given by

Λo2p2p Z2p2p  Z2p2p Λc2p2p ¼  Λo2p2p XTon2p Bnm Γm2p

ð4Þ

Note that the Sylvester equation presented in Eq. (4) implements the QPEVA technique. In Eqs. (3) and (4), 2p Γ ¼ fγ i g2p i ¼ 1 A Cm2p is known as parametric matrix. The vector set fγ i gi ¼ 1 is closed under complex conjugation in the same 2p order as the set μi i ¼ 1 . The deterministic MNQPEVA problem is then given by min Γ

 1 ‖FðΓÞ‖2F þ ‖GðΓÞ‖2F 2

ð5Þ

where J : J F denotes the Frobenius norm of a matrix. 2.2. Brief comparison of MNQPEVA with H 1 A comparison of the deterministic MNQPEVA technique (Eqs. (3)–(5)) with the well established H 1 technique [17,25,26,37–49] is presented in this section using a 2 and a 500 degree of freedom systems. The comparison is based on the computation time and the actuator capacity required by the two techniques for design/implementation of closed-loop

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

363

systems for the same open-loop system (viz. the 2 and 500 degree of freedom systems). To compare the computation time, all computations in this section are performed in 32-bit MATLAB [61] using a single core in a computer with Intel(R) Core (TM)2 Duo processor having speed of 2 GHz and RAM of 3 GB. The mathematical formulations associated with the H 1 technique used in this section can be found in the literature [5, Section 10.6] and not presented here to save space. The 2 degree of freedom system, considered for comparison of the two techniques, is shown in Fig. 3. The masses m1;2 ¼ 104 kg and the stiffness k1 ¼ 0:7 MN=m and k2 ¼ 1:2 MN=m. Rayleigh type damping matrix with ξ1;2 ¼ 0:02, where ξ1;2 denote the damping ratios associated with the first two modes of the system, is assumed for this example problem. The eigenvalues of the open-loop system (viz. the 2 degree of freedom system) are presented in Table 1. The system is subjected to a scaled ground motion g€ which is recorded at the Rinaldi receiving station (fault-normal component) during the 1994 Northridge earthquake with moment magnitude, Mw ¼6.7. The open-loop system exhibits high dynamic response (see Fig. 4) as the dominant range of the loading spectrum is around the first damped natural frequency of the system (5.476 rad/ s) associated with the first pair of eigenvalues (see Table 1). The dynamic response of the 2 degree of freedom system can be reduced using both the deterministic MNQPEVA and H 1 control techniques. For this purpose, the control matrices for both the techniques are chosen in such a manner that they are equivalent and comply with the location of an actuator at floor level 1 (see Fig. 3). It can be observed from Table 1 that the H 1 technique shifts all the eigenvalues of the open-loop system. Since the MNQPEVA technique is capable of shifting a chosen set of open-loop eigenvalues to a target set, two cases are considered here while designing a closed-loop system using the MNQPEVA technique – (i) all eigenvalues of the open-loop system are shifted and matched to the eigenvalues of the H 1 based closed-loop system and (ii) only the first pair of eigenvalues associated with the first damped natural frequency of the open-loop system are shifted. It can be observed from Fig. 4 and Table 1 that the closed-loop systems designed

Fig. 3. Schematic representation of a two degree of freedom system subjected to ground motion g€ . Table 1 Comparison of deterministic MNQPEVA and H1 techniques for a 2 degree of freedom system. System

Eigenvalues

Computation time (s)

Open-loop H1 closed-loop MNQPEVA closed-loop matched to H1 MNQPEVA closed-loop

f  0:110 7 5:476i;  0:335 7 16:730ig – f  0:484 7 5:464i;  1:273 7 16:755ig 2.086 f  0:484 7 5:464i;  1:273 7 16:755ig 0.801 f  1:550 7 6:003i;  0:335 7 16:730ig 0.149

Optimization algorithm

Actuator force capacity (kN)

J ΔK J F JKJF

– Interior-point Quasi-Newton line search Quasi-Newton line search

20.43 20.43

1.73 1.73

376.75 376.75

59.57

5.00

681.77

(%)

J ΔC J F JCJF

(%)

Note: ΔK and ΔC denote the changes in the stiffness and damping matrices, respectively, of the closed-loop system with respect to their open-loop counterparts.

Fig. 4. Displacement comparison at floor level 1 for the open-loop system and closed-loop systems designed according to the deterministic MNQPEVA and H1 techniques.

364

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

using the H 1 technique and the MNQPEVA technique matched to H 1 based closed-loop system yield same results in terms of the response and the required actuator capacity, however, they differ in computation time. The closed-loop system designed using the MNQPEVA technique with only the first pair of open-loop eigenvalues shifted illustrates the advantage of the MNQPEVA technique over the H 1 technique in terms of achieving better control on the response of the system with much lesser computation time (see Fig. 4 and Table 1). It should, however, be noted that the better closed-loop system is achieved at the expense of higher actuator capacity (see Table 1). The effect of the size of an open-loop system on the computation time incurred by the deterministic MNQPEVA and H 1 techniques is also investigated and presented here using a 500 degrees of freedom large-scale system with m1;…;500 ¼ 104 kg, k1;…;500 ¼ 25 GN=m, and Rayleigh type damping matrix with ξ1;2 ¼ 0:02. As expected, the H 1 technique shifts all the openloop eigenvalues. For this example problem, the MNQPEVA technique shifts only the first pair of open-loop eigenvalues associated with the first damped natural frequency of the system that lies in the dominant range of the loading spectrum. For this problem, the interior-point and quasi-Newton line search optimization algorithms are used to implement the deterministic H 1 and MNQPEVA techniques, respectively. It is observed that the computation time required by the deterministic MNQPEVA and H 1 techniques is 52.35 s and 899.78 s, respectively, rendering the MNQPEVA technique to be more suitable for the design of controllers for the large-scale systems. It should be remarked here that the difference in computational time in the two examples presented above depends on several factors, for example, formulations of the optimization problems (e.g., constrained H 1 problems, unconstrained MNQPEVA problems), size of the system defined by the degrees of freedom, and the underlying space on which the problem needs to be constructed (i.e., state space based H 1 problem, physical space based MNQPEVA problem). The results presented in this section illustrate the drawbacks associated with the H 1 technique mentioned earlier (see Section 1) and highlight the need for the development of a computationally effective stochastic QPEVA technique for economic and resilient design of the closed-loop systems.

3. Formulation of methodology for stochastic closed-loop systems In this section, a variant of the deterministic MNQPEVA problem which minimizes the norm of feedback matrices as well as the mass of the closed-loop system is presented. This variant problem is termed as minimum norm quadratic partial eigenvalue assignment-minimum mass closed-loop system (MNQPEVA-MMCLS) problem. The mass of the closed-loop system is minimized in order to reduce the material cost and the inertial load induced by dynamic effects. Finally, based on the variant of the deterministic MNQPEVA problem, an OUU is formulated for economic and resilient design of stochastic closed-loop systems. In the deterministic MNQPEVA-MMCLS problem, structural (material and geometric) and control input parameters are considered as decision variables. Also, the response quantities (i.e., displacement, velocity, acceleration, stress, etc.) are constrained within certain limit to ensure safety. The deterministic MNQPEVA-MMCLS problem is described as follows: !2  g ðqÞ MðqÞ 2 min a  1 þ ð1  aÞ q g target Mavg s:t:

yi ðqÞ ryicr ; i ¼ 1; …; ny

lqj rqj r uqj ; j ¼ 1; …; nq

ð6Þ

where s.t. refers to “subject to”, q is a vector of structural and control parameters, a is a positive real-valued weight coefficient, MðqÞ is the total mass of the closed-loop system for given q, and gðqÞ is given by  h     i 1 J F q; Γ J 2F þ J G q; Γ J 2F ð7Þ g ðqÞ ¼ min 2 Γ which is same as the formulation of MNQPEVA problem [32–34]. The terms gtarget and Mavg in Eq. (6) are normalizing factors used to eliminate any bias due to difference in magnitude of the two terms while numerically carrying out the optimization scheme. In addition, gtarget serves as a target value for g ðqÞ. In Eq. (6), yi ðqÞ represents the i-th response quantity (e.g., displacement, velocity, acceleration, and stress) with the associated critical value given by yicr . The terms lqj and uqj represent the lower and upper bounds, respectively, on the j-th decision variable qj. The terms ny and nq denote the number of response quantities and decision variables, respectively. Formulation of stochastic MNQPEVA-MMCLS problem assumes a random vector q which collectively represent all random parameters associated with the structural system and control input. Using the concept of OUU [56], the stochastic MNQPEVA-MMCLS problem can be written as 2 3 !2  g ðqÞ MðqÞ 2 5 min E4a 1 þ ð1 aÞ μq g target Mavg s:t:

P½yi ðqÞ r yicr  ¼ 1  P f y ZP yi ; i ¼ 1; …; ny i

P½lqj r qj ruqj  ZP qj ; j ¼ 1; …; nq

ð8Þ

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

365

Fig. 5. Pictorial presentation of constraint satisfaction for probability of failure: (a) existing methods and (b) proposed method.

where P yi and P qj are the user-defined lowest probabilities to be satisfied for the i-th constraint and the j-th input bounds, respectively. It should be noted that ð1  P yi Þ is the prescribed probability of failure for yi ðqÞ, where failure is defined by yi ðqÞ 4yicr . The term P f y denotes the current probability of failure for yi ðqÞ given the vector of design variables μq which i represents the mean vector of q for the MNQPEVA-MMCLS problem defined above by Eq. (8). The solution of an OUU problem is typically obtained by solving a simplified version, known as robust optimization (RO) problem, where the objective function and the constraints are expressed as weighted sum of mean and variance of respective quantities, neglecting the higher order statistics [62–66,56]. The available algorithms for solving RO are restrictive because (i) they typically assume statistical independence among the random input parameters and (ii) they require a priori knowledge of the probability density functions (pdfs) of constraints and objective function for proper selection of weight coefficients for mean and variance. In RO, the probabilistic constraint associated with the i-th random response yi ðqÞ is satisfied by shifting the support of the pdf pyi ðqÞ of yi ðqÞ to the left of yicr without changing the shape of pyi ðqÞ , where yicr represents the critical value of yi ðqÞ indicating failure (see Fig. 5(a)). The constraint on P f y can, however, be satisfied via a i combination of shift in support and change in shape of pyi ðqÞ (see Fig. 5(b)) by considering the set of distribution parameters nθ fθk gk ¼ 1 of the joint pdf pq of q as the set of design variables for the MNQPEVA-MMCLS problem. The term nθ denotes the number of distribution parameters required to define the joint pdf pq . This idea can be implemented by using a family of pdfs that enjoy flexibility in shape (e.g., Beta distribution). It should, however, be noted that the proposed idea uses only the distribution parameters and is in general independent of the family of distributions. A designer should use engineering judgement to select the appropriate family of pdfs complying with the manufacturing/construction specifications. Using the proposed idea, the stochastic MNQPEVAP-MMCLS problem can be formulated as follows: 2 !2    !2 3 nθ nθ g qjf θ g M qjf θ g k k k¼1 k¼1 5 E 4a  1 þ ð1  aÞ min n g target Mavg fθk g θ k ¼ 1

s:t:

P½yi ðqjfθk gk θ¼ 1 Þ ryicr  ¼ 1 P f y Z P yi ; i ¼ 1; …; ny n

i

lθk r θk r uθk ; k ¼ 1; …; nθ ð9Þ     nθ nθ where g qjfθk gk ¼ 1 , M qjfθk gk ¼ 1 are the same as mentioned earlier in accordance with Eq. (6), however, now expressed n as functions of fθk gk θ¼ 1 . The terms lθk and uθk represent the lower and upper bounds on the k-th decision variable, respectively.   n The computational bottleneck of this new formulation (Eq. (9)) is efficient yet accurate estimation of P f y fθk gk θ¼ 1 ; 8 i.   i nθ This is due to the reason that the accuracy of the estimate of P f y fθk gk ¼ 1 ; 8 i depends on the knowledge of failure domain i Sif ¼ yi ðqÞ Z yicr which might be disjointed and complicated in shape, and not known explicitly as a function of system parameters (for example, available only through large-scale finite element simulations). To mitigate this computational   n difficulty associated with accurate and efficient estimation of P f y fθk gk θ¼ 1 ; 8 i; a novel scheme based on hierarchical i clustering technique (HCT) and importance sampling (IS) scheme has been developed and presented in the following section.

4. Estimation of probability of failure The optimality of the solution of OUU problem depends on the accuracy of the estimated probability of failure Pf. Estimation of Pf requires the knowledge of failure domain Sf. Accuracy and efficiency in estimation of Pf is compromised when Sf is disjointed and complicated in shape (see Fig. 6), and not known explicitly as a function of system parameters. The first step to increase the accuracy and efficiency of the estimate of Pf is to reliably identify the failure domain Sf. In this regard, estimation of Pf involves two tasks: (i) identification of disjointed and complicated shaped Sf and (ii) computation of Pf based on the samples digitally simulated form the identified Sf. These two tasks are discussed below in detail.

366

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Fig. 6. Schematic representation of (a) variation of yðqÞ with q and (b) disconnected Sf associated with yðqÞ where q ¼ ½x1 ; x2 .

4.1. Identification of failure domain Identification of the unknown Sf can be initiated by applying clustering techniques on available digitally simulated failure samples even if they are only few in numbers. Towards this identification process, the hierarchical clustering technique (HCT) is preferred over the k-means and other clustering techniques in this work because HCT does not assume the failure samples to be independent and identically distributed [60, Section 14.3.12]. The HCT is based on two approaches:

 agglomerative and divisive (AHCT and DHCT). Time complexities associated with the AHCT and DHCT are O N 3 and O 2N , respectively, where N is the number of samples. This implies that the AHCT is more efficient compared to the DHCT and should be preferred. In AHCT, the samples are sequentially grouped together based on (i) a distance measure among the samples (e.g., Euclidean, Minkowski, and Cosine) and (ii) a distance measure among clusters known as linkage (e.g., Single linkage (SL), Complete linkage (CL), and Average linkage (AL)) [60, Section 14.3.12]. The process followed in AHCT is, in
 general, pictorially depicted via a dendrogram. Note that the time complexity of the SL based AHCT is O N2 whereas that
 of the CL and AL based AHCTs are O N 2 logðNÞ [67]. Recent literatures also show that the SL based AHCT is more stable and has better convergence criteria [68]. Therefore, the SL based AHCT is employed in this work to identify the unknown Sf. The algorithmic steps are presented in Algorithm 1. Algorithm 1. Identification of disjointed failure domain. Input: n-D Euclidean domain representing support S of q, failure criteria defined by yicr , and exploratory sample budget Nex for domain exploration. Step 1: Obtain Nex exploratory samples in S assuming elements of qj to be statistically independent and qj  U½Sminqj ; Smaxqj , 8 j to explore the failure domain Sf A S without any bias. Note that Sminqj and Smaxqj are the possible minimum and maximum bounds on qj , respectively, to be chosen for numerical purpose. Step 2: Segregate the failure samples to identify the failure domain Sf. Step 3: Apply the SL based AHCT to obtain a dendrogram containing information about the disjointedness of the Sf. Step 4: Identify nSf parts of a disjointed Sf by slicing the dendrogram using a cut-off value associated with the normalized incremental height of the dendrogram. Output: Disjointed parts of Sf.

To illustrate the proposed idea, an example problem with 2 random variables x1 and x2 is presented here. For this problem, random vector q ¼ ½x1 ; x2 T and the support S of q is considered as S ¼ ½0; 1  ½0; 1. The joint pdf of q is assumed to be pq ðqÞ ¼ px1 ðx1 Þpx2 ðx2 Þ with x1  Betað4; 3; 0; 1Þ and x2  Betað6; 2; 0; 1Þ. The joint pdf is depicted in Fig. 7(a). The failure 2 2 2 2 domain Sf comprises two disconnected parts, Sf 1 ¼ fq: ðx1  0:090Þ þ ðx2  0:120Þ r1g and Sf 2 ¼ fq: ðx1  0:900Þ þ ðx2  0:070Þ r1g (see 0:0752 0:1002 0:0752 0:0502 Fig. 7(b)). The true probability of failure for this example problem is 3:763  10  6 . Following Step 1, N ex ¼ 500 exploratory samples in S have been generated to explore the disjointed Sf (see Fig. 7(b)). This results in N f ¼ 16 failure samples. The failure samples are segregated (Step 2) and applying the SL based AHCT two clusters of failure samples are obtained using a cut-off value of 0.1 for the normalized incremental height of the dendrogram (Steps 3–4) (see Fig. 7(c)). The clusters of samples well identify the true failure domain (see Fig. 7(d)). 4.2. Computation of probability of failure Once Sf is identified, the conventional MCS technique can be employed to estimate the Pf associated with a response quantity y [57,58]. The MCS technique, however, requires sufficiently large number of samples for accurate estimation of low Pf typically associated with the tail region of pdf of y, hence, computationally inefficient. To mitigate this difficulty, IS scheme is typically employed since it requires less samples. An IS estimate of an order statistic of f ðxÞ representing a function of the random variable x is obtained by first generating samples from the domain of interest SI (i.e., Sf in the current context) using an auxiliary pdf g x ðxÞ defined over SI and then correcting for the bias using original pdf px ðxÞ [57–59]. The choice of g x ðxÞ is governed by (i) 8 x A SI , g x ðxÞ 4 0 if px ðxÞ 4 0 and (ii) px ðxÞ okg x ðxÞ with k⪢1 to ensure almost sure convergence of the IS scheme [58]. In addition, if R 2 p2x ðxÞ SI f ðxÞgx ðxÞ dx o 1 holds then the estimator has finite variance [58]. At this stage one should note that the traditional IS scheme cannot be used to estimate Pf if information about “exact” Sf is unavailable. To relax this restriction, a modified version of the

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

367

Fig. 7. Various attributes of the example problem: (a) probability density function, (b) failure domain, (c) sliced dendrogram, and (d) clustered failure samples.

Fig. 8. Actual and importance distributions: (a) traditional IS scheme and (b) modified IS scheme.

traditional IS scheme, where samples are generated from a domain SS + SI (or Sf), is proposed in this paper. Let x be a random variable and we are interested in the estimation of P½x A Sf . The IS estimate of P½x A Sf  is given by N SS   1 X px ðxðiÞ Þ IS P x A Sf ¼ NSS i ¼ 1 g x ðxðiÞ Þ f

ð10Þ

where N SS is the number of samples generated within SS. To illustrate the efficiency of this IS scheme (Eq. (10)) over the MCS techniques, x  Rayleighð1Þ with S ¼ ½0; 1Þ is con2 sidered. The pdf of x is given by px ðxÞ ¼ xe  x =2 IS where IS is the associated indicator function (see Fig. 8(a)). Let Sf be ½3:2; 4:0. Analytically, P½x A Sf  is 0.00564. The MCS estimate of P½x A Sf  with 105 samples is 0.00596 with 5.66% error. Using the IS scheme given in Eq. (10) with 500 samples, the estimate of P½x A Sf  is 0.00556 with 1.35% error. It should, however, be noted that the auxiliary pdf g x ðxÞ ¼ ð4:5 1 2:7ÞISS with SS ¼ ½2:7; 4:5 has been used for this IS scheme (see Fig. 8). When the failure domain Sf consists of disconnected parts with their shapes approximately known, the estimate of   probability
of  failure P½qA Sf  can be obtained by applying the modified IS scheme repeatedly for each disconnected failure domain Sf j and summed over all j. The algorithmic steps associated with such estimate are presented in Algorithm 2.

368

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Algorithm 2. Computation of probability of failure. Input: Identified failure domain Sf, tolerance ϵ, joint pdf pq ðqÞ, and sample budget N. Step 1: Obtain n-D boxes bounding each of the disjointed part of Sf with tolerance ϵ. Step 2: Obtain a total of N new samples from the bounding boxes using auxiliary pdfs gjq ðqÞ; 8 j. The sample size Nj associated with each box is proportional to the size of the box.  p ðqðiÞ Þ PnS PN Step 3: Use the IS scheme P f ¼ j ¼f 1 N1j i ¼j 1 f qðiÞ qj ðiÞ ISf to estimate probability of failure. gq ðq Þ

j

Output: IS estimate of probability of failure Pf.

Applying Algorithm 2 on the 2D example problem described in Section 4.1, the IS estimate of Pf is obtained as 3:696  10  6 with 1.77% error for N ¼600 samples and ϵ ¼ 0:1. The boxes bounding the disjointed parts of failure domain (Step 1) as well as the N ¼600 samples generated within these boxes (Step 2) are depicted in Fig. 9. To assess the accuracy of the proposed scheme, a probabilistic convergence study is conducted with 1000 different sets of N ex ¼ 500 exploratory samples (Step 1 of Algorithm 1). The convergence study indicates that the 95% quantile range, mean, and standard deviation of the absolute error in Pf reduce significantly with sample budget N (see Fig. 10(a, b)). The convergence study also reveals that the 95% quantile range and the mean of the computation time increases almost linearly with sample budget N (Fig. 11).

5. Recommendation of design guidelines  n

The methodology developed in Section 3 results in optimal values of distribution parameters fθk gk θ¼ 1 for the structural and control input properties. Probabilistic characterization of the elements of feedback matrices F and G associated with the optimal

Fig. 9. Samples generated from boxes bounding the failure clusters (2D example problem).

Fig. 10. Convergence of the novel scheme used to estimate Pf (2D example problem).

Fig. 11. Computation time required by the novel scheme to estimate Pf (2D example problem).

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

369

Fig. 12. Plots depicting (a) the 2D steel frame structure and (b) a few FRFs of the frame along with the normalized smooth FA of the ground motion.

Table 2 Lower and upper bounds of distribution parameters associated with pdfs of ρ and E. ρ

Parameter

E

Distribution parameter

α1

β1

cl1 (kg m

Lower bound Upper bound

0 1

0 1

7000 7320

3

)

cu1 (kg m 8280 8600

3

)

α2

β2

cl2 (GPa)

cu2 (GPa)

0 1

0 1

180 192

228 240

 n

stochastic closed-loop system can be obtained using fθk gk θ¼ 1 . However, for implementation purpose, a unique set of feedback matrices is required and should be selected from the support of probabilistically characterized set of feedback matrices. Such a selection criterion for the unique set of feedback matrices reasonably ensures the resilience level selected for the design of the closed-loop system. Additional conditions can be imposed on the selection process for the unique set of feedback matrices so that the chosen set guarantees robustness in the response of the designed system. Towards this end, an optimization problem has been formulated in this work that results in an optimal unique set of feedback matrices which guarantees the selected resilience level as well as the robustness in response of the designed system. The optimization problem is given by h i min E ðyðM; C; K; B; F; GÞ  E½yðM; C; K; B; F; GÞÞ2 F;G

s:t:

lF j;k r F j;k r uF j;k ; lGj;k rGj;k r uGj;k ; j ¼ 1; …; m; k ¼ 1; …; n

ð11Þ

where yðM; C; K; B; F; GÞ is the response with respect to which robustness is sought. The terms lF j;k , uF j;k , and lGj;k , uGj;k denote the bounds on the (j,k)-th element of the matrices F and G, respectively. It should be remarked here that the bounds are associated with the support of the marginal pdfs of the elements of the matrices F and G. This optimization problem will involve a large number of decision variables because the number of elements in F and G will be some integer multiple of the number of nodes in the FE model of the large-scale structure. Such a problem can be efficiently solved using evolutionary algorithms [69, Chapter 4].

6. Numerical illustration The methodology developed in this paper (Sections 3-5) for design of stochastic closed-loop system is presented here numerically using a 4-story steel building (similar to those considered in the literature [70,71]; see Fig. 12(a)) subjected to a ground acceleration recorded at the Rinaldi receiving station (fault-normal component) during the 1994 Northridge earthquake with moment magnitude, Mw ¼6.7. The acceleration record is scaled to obtain a peak ground acceleration (PGA) of 0.4 g where g is the acceleration due to gravity. In this problem, only the structural material properties (density ρ and Young's modulus E) are assumed to be random while keeping the dead-load (w), the damping ratios associated with the first two modes of vibration (ξ1 and ξ2), and all the 3 3 3 geometric properties (fbi gi ¼ 1 ; fhi gi ¼ 1 ; ft f i g3i ¼ 1 ; ft wi g3i ¼ 1 , and fli gi ¼ 1 ) as deterministic. Note that the terms bi ; hi ; t f i ; t wi , and li represent the breadth, height, flange thickness, web thickness, and length of the I-sections used as structural elements of the frame (see Fig. 12(a)), respectively. One can, however, assume all the parameters to be random at the cost of increase in dimensionality. The values of the deterministic parameters assumed in this example are w ¼ 2800 kg=m, ξ1 ¼ ξ2 ¼ 0:02, b1 ¼ 0:7754 m, h1 ¼ 0:3190 m, t f 1 ¼ 0:0229 m, t w1 ¼ 0:0381 m, b2 ¼ 0:4922 m, h2 ¼ 0:3781 m, t f 2 ¼ 0:0330 m, t w2 ¼ 0:0445 m, b3 ¼ 0:4826 m, h3 ¼ 0:3938 m, t f 3 ¼ 0:0254 m, t w3 ¼ 0:0381 m, l1 ¼ 7:31 m, l2 ¼ 5:49 m, and l3 ¼ 3:66 m. It is assumed that     ρ  Beta α1 ; β1 ; cl1 ; cu1 , E  Beta α2 ; β2 ; cl2 ; cu2 , and ρ, E are statistically independent. The lower and upper bounds of the
 2 are tabulated in Table 2. Assuming α1 ¼ 3, β1 ¼ 2, distribution parameters fαj g2j ¼ 1 ; fβj gj ¼ 1 ; fclj g2j ¼ 1 ; fcuj g2j ¼ 1

370

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Table 3 Lower and upper bounds of distribution parameters associated with pdfs of r and m. Parameter

r

m

Distribution parameter

α3

β3

cl3

cu3

α4

β4

cl4 (rad s  1)

cu4 (rad s  1)

Lower bound Upper bound

0 1

0 1

 1.00  0.88

 0.52  0.40

0 1

0 1

16.0 19.2

28.8 32.0

Fig. 13. Exploration of disjointed failure domain: (a) samples for exploration (blue ¼safe and red¼failure) and (b) the disjointed parts. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Table 4 Representative initial and optimal values of distribution parameters associated with pdf of ρ. Distribution parameter

α1

β1

cl1 (kg m  3)

cu1 (kg m  3)

Representative initial Optimal

3 0.0066

2 5.7122

7000 7001.1

8600 8390.4

cl1 ¼ 7000 kg=m3 , cu1 ¼ 8600 kg=m3 , α2 ¼ 2, β2 ¼ 3, cl2 ¼ 180 GPa, cu2 ¼ 240 GPa, and the deterministic values mentioned above as the initial design of the stochastic structure (in this case the 2D frame), it is found that the fundamental damped natural frequency ωdn1 of the structure lies within the dominant range of the spectrum of the ground motion (see Fig. 12(b)). For AVC, fλ1;2 g associated with ωdn1 of the structure need to be replaced by new target eigenvalues fμ1;2 ¼ r 7 jmg where pffiffiffiffiffiffiffiffi j ¼ 1. To account for uncertainty in control input, fr ¼ Reðμi Þ; m ¼ Imðμi Þ; i ¼ 1; 2g have been assumed to be random with     r
 Beta α3 ; β3 ; cl3 ; cu3 and m Beta α4 ; β4 ; cl4 ; cu4 . The lower and upper bounds of the distribution parameters 4 4 4 4 fαj gj ¼ 3 ; fβj gj ¼ 3 ; fclj gj ¼ 3 ; fcuj gj ¼ 3 are tabulated in Table 3.

For this example problem, the random vector q ¼ ½ρ; E; r; mT with ρ, E, r, m o i.e., n being statistically independent, 16 4 pq ðqÞ ¼ pρ ðρÞpE ðEÞpr ðrÞpm ðmÞ, and the set of design variables is given by fθk gk ¼ 1 ¼ fαj g4j ¼ 1 ; fβ j gj ¼ 1 ; fclj g4j ¼ 1 ; fcuj g4j ¼ 1 . The control matrix B for this problem is chosen as " #T B1;1 ¼ 1 0 … 0 0 0 … 0 B1002 ¼ 0 0 … 0 B2;52 ¼ 1 0 … 0

to comply with the placement of actuators at the locations associated with Node-1 and Node-52 (see Fig. 12(a)). For design, the maximum von-Mises stress, denoted by σ VMmax , is considered as the only response quantity y1 of interest. Failure occurs when y1 exceeds the yield stress of steel denoted as σyield. It is assumed that σ yield ¼ 250 MPa. The maximum possible support of q ¼ ½ρ; E; r; mT obtained using the lower bound values of clj 4j ¼ 1 and the upper bound values of cuj 4j ¼ 1 has been explored using Nex ¼ 50 000 samples and three disjointed parts of failure domain are detected (Algorithm 1). The disjointed failure domain is depicted in Fig. 13 in the E–r–m space. Note that the sample budget N ¼ 50 000 has been considered in this example for computation of Pf when q assumes the maximum possible support (Algorithm 2). To solve the stochastic MNQPEVA-MMCLS problem given in Eq. (9), a resilience level associated with P y1 ¼ 0:9999 is selected. The coefficient a is chosen as 0.5 to provide same weight to both the terms. For this example, the values of the normalization factors gtarget and Mavg are assumed to be 8:3107  1013 and 42 829, respectively. The optimization problem is solved in MATLAB [61] using genetic algorithm (GA) followed by sequential quadratic programming (SQP). The GA is employed for this problem with 0.8 crossover fraction, 0.01 mutation rate, 1 elite count, 25 initial population, and tour4 nament selection. A representative initial and the optimal values of fθk gk ¼ 1 ¼ fα1 ; β 1 ; cl1 ; cu1 g (associated with ρ),

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

371

Table 5 Representative initial and optimal values of distribution parameters associated with pdf of E. Distribution parameter

α2

β2

cl2 (GPa)

cu2 (GPa)

Representative initial Optimal

2 3.2673

3 5.3946

180 180

240 240

Table 6 Representative initial and optimal values of distribution parameters associated with pdf of r. Distribution parameter

α3

β3

cl3

cu3

Representative initial Optimal

4 3.6125

3 6.1900

 1.00  0.8866

 0.40  0.4317

Table 7 Representative initial and optimal values of distribution parameters associated with pdf of m. Distribution parameter

α4

β4

cl4 (rad s  1)

cu4 (rad s  1)

Representative initial Optimal

1 6.4085

1 0.0629

16.0 17.3873

32.0 30.6424

Fig. 14. Normalized histograms of (a) representative initial and optimal closed-loop systems and (b) optimal and sub-optimal closed-loop systems.

Table 8 Estimate of various statistics of σ VMmax for various stochastic systems obtained through solving various optimization problems. Case

Mean

Standard deviation

Coeff. of variation

Initial closed-loop Optimal closed-loop Sub-optimal closed-loop Optimal feedback Unique feedback

172.83 166.52 162.81 121.03 164.75

34.62 14.03 29.32 2.59 4.55

0.2003 0.0843 0.1801 0.0214 0.0276

fθk gk ¼ 5 ¼ fα2 ; β2 ; cl2 ; cu2 g (associated with E), fθk gk ¼ 9 ¼ fα3 ; β3 ; cl3 ; cu3 g (associated with r) and fθk gk ¼ 13 ¼ fα4 ; β 4 ; cl4 ; cu4 g (associated with m) are presented in Tables 4, 5, 6 and 7, respectively. Note that the probability of failure for the optimal closed-loop system P opt ¼ 0:0000 as compared to the probability of failure for the initial closed-loop system P ini f y ¼ 0:0247. fy 1 1 The distributions of σ VMmax of the closed-loop systems, in the form of normalized histograms obtained from 8

12

16

10000 samples generated using the values of fθk gk ¼ 1 before and after optimization, are presented in Fig. 14(a). It can be observed from this figure that the representative initial closed-loop systems cannot withstand the dynamic load with required level of safety in the presence of uncertainties, whereas, their counterpart closed-loop systems can. Estimates of various statistical quantities associated with σ VMmax for these different cases are presented in Table 8. It can be noted from Table 8 that the standard deviation of σ VMmax for the optimal closed-loop system is quiet less and hence good from design view-point. At this stage, a designer may be concerned about how to maintain the optimal probability distributions 16

372

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

of ρ and E during the fabrication/construction of structural members since the state-of-the-art manufacturing process and/ or construction technology are limited in their capabilities to produce material and geometric parameters with target distributions. Also, maintaining the optimal probability distributions of r and m might also be difficult. To get some insight on how the designed stochastic closed-loop system will behave under such circumstances, qj  Uðcopt ; copt uj Þ; 8 j have been lj opt opt opt opt opt opt assumed (i.e., ρ  Uðcopt ; copt u1 Þ; E  Uðcl2 ; cu2 Þ; r  Uðcl3 ; cu3 Þ, and m  Uðcl4 ; cu4 Þ) and another 10000 samples of q are l1

generated. Note that the terms copt and copt uj denote the optimal lower and upper bounds of the support of the random lj variable qj , respectively. The choice of uniform distribution with the optimal bounds for the random parameters comes from the fact that a uniform distribution has maximum entropy and can account for the maximum variability within these bounds. A comparison of the normalized histograms of the optimal and the new sub-optimal closed-loop systems is depicted in Fig. 14(b). One can observe from this figure that the mean of σ VMmax for the optimal and sub-optimal cases is almost same while the standard deviation of σ VMmax for the sub-optimal case is higher compared to the optimal counterpart. ¼ 0:0104 which exceeds the user prescribed It is observed that the probability of failure for the sub-optimal case P subopt f y1

value. The sub-optimal case, nevertheless, still represents a significant improvement over the representative initial closedloop system in restricting the response (see Fig. 14(b)). This is practically very appealing since the proposed methodology can be confidently adapted in the design of active controlled system to account for the effects of uncertainties. To investigate the effect of weight coefficient a on optimal results, the optimization problem given in Eq. (9) is solved again with a¼1. Note that the objective function of this OUU problem does not have any contribution from MðqÞ. Due to this reason, the solution for this case will provide the optimal distribution of feedback matrices only and not the optimal distribution of structural properties (material and geometric). This approach (i.e., solution of Eq. (9) with a¼1.0) will be practically appealing since the state-of-the-art manufacturing process is currently limited in their capability to produce material and geometric parameters, and structural complexities with target optimal pdfs. This choice of the weight coefficient ðaÞ will allow the current manufacturing technologies to embrace the design methodology proposed in this paper. The results for the a ¼1 case are presented in Fig. 15 and Table 8. It can be observed from Table 8 that the mean and standard deviation of σ VMmax for a¼1 case are much lesser (about 27 % for mean and 81% for standard deviation) compared to that of a ¼0.5 case. However, for both a ¼0.5 and a ¼1 cases, the mean values of gðqÞ are very close to the gtarget with percent differences in the order of 10  4.  n After the optimal distribution parameters fθk gk θ¼ 1 of q are obtained, one can find out the distributions of the elements of stochastic feedback matrices F and G by generating samples of q using MCS technique first and then applying the deterministic MNQPEVA technique to find the feedback matrices F and G for all the samples. Utilizing a set of 10 000 samples  16 generated using the optimal distribution fθk gk ¼ 1 of the stochastic closed-loop system parameters and fitting Beta distributions to the resulting samples of F and G, estimates of the distribution parameters of the elements of F and G have been obtained (results not reported). A high correlation among the elements of F and G is observed, which is expected. Using these information, one can solve the optimization problem given in Eq. (11) to obtain a unique set of feedback matrices for implementation purposes. Evaluation of the objective function for this optimization problem will, however, require large number of samples of system matrices M; C; K (conventional MCS technique) thereby rendering the process to be computationally inefficient. Due to this reason, sub-optimal solution is obtained by (i) considering 500 samples of M; C; K and 500 samples of F; G obtained from the 10 000 samples mentioned above using resampling technique and (ii) selecting the set of F; G among the 500 samples that led to the minimum value of the objective function given in Eq. (11). The results for this chosen F and G have been presented in Fig. 16 and Table 8. The results indicate that the mean of σ VMmax for this case is similar to that obtained for the optimal closed-loop system case (with a¼0.5). Also, the probability of failure associated with the unique set of feedback matrices P feedback ¼ 0:0000 rendering it very appealing for implementation purposes. fy 1 An indirect comparison of the cost associated with the various design cases presented in Table 8 for active control of stochastic systems is also shown here in terms of the 95% quantile range for the actuator capacity (Fig. 17) and the 95% quantile range for the total mass of the stochastic closed-loop system (Fig. 18). Note that the actuator capacity has been represented in terms of the peak control force, peak velocity, and peak displacement of the actuator. It can be observed from Figs. 17 and 18 that the optimal closed-loop system, sub-optimal closed-loop system, closed-loop system with optimal feedback, and closed-loop system with unique feedback require actuators with less capacity and involve less mass as compared to the representative initial closed-loop system. One can also observe from Fig. 17 that the peak control force associated with the closed-loop system with optimal feedback matrices is much less compared to the peak control force required for the optimal closed-loop system. The trend for the total mass associated the closed-loop system with optimal feedback matrices and the optimal closed-loop system is, however, reversed (see Fig. 18). Such trends comply with the fact that for the optimal closed-loop system both the feedback gain and the mass are minimized whereas for the closed-loop system with optimal feedback matrices only the feedback gain is minimized. It can be also seen from Figs. 17 and 18 that the 95% quantile range for the actuator capacity and the total mass associated with the closed-loop system with unique feedback matrices are always contained within the 95% quantile range for the respective quantities associated with the optimal closed-loop system. This is because the unique feedback matrices ðF; GÞ are obtained based on the samples of F and G that  16 are generated using the optimal distribution fθk gk ¼ 1 associated with the optimal closed-loop system parameters.

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

373

Fig. 16. Normalized histograms of representative initial closed-loop system, optimal closed-loop system, and closed-loop system with unique feedback matrices.

Fig. 17. 95% quantile range for actuator capacities in terms of (a)–(b) peak control force, (c)–(d) peak velocity, and (e)–(f) peak displacement of the actuators located at Node-1 and Node-52, respectively, for the different cases presented in Table 8.

7. Conclusions A new methodology for economic and resilient design of active controlled large-scale stochastic structural systems has been formulated in the present work. The methodology is based on the concepts of QPEVA technique and optimization under uncertainty. Accurate yet efficient estimation of low probability of failure associated with disjointed and complicated shaped failure domain is the computational bottleneck for the proposed methodology. A novel scheme based on hierarchical clustering technique and importance sampling scheme has been proposed here to mitigate the computational difficulty associated with estimation of this probability of failure. Following the methodology, a designer can obtain the optimal

374

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

Fig. 18. 95% quantile range for total mass of the system for the different cases presented in Table 8.

probability distributions of the random input parameters (structural and control). Concerns that might arise from fabrication/construction view-point about maintaining the optimal probability distributions of the random input parameters have also been discussed and two way-outs in the form of (i) sub-optimal closed-loop systems and (ii) closed-loop systems with optimal feedback matrices have been presented. The second way-out is slightly conservative as compared to the first one but better from design view-point as (i) it ensures required resilience level and (ii) it leads to more economic design in terms of the energy required for control purposes. A technique to obtain a unique set of feedback matrices from the stochastic feedback matrices is also recommended in this paper for implementation purposes. The resulting unique set of feedback matrices ensures the required resilience level and hence very appealing from design perspective.

Acknowledgements The authors acknowledge the support of computational resources provided by the Center for Computational Research at the University at Buffalo (UB) and partial financial support under the UB RENEW seed grant 1120795-1-64755.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

K. Billah, R. Scanlan, Resonance, tacoma narrows bridge failure, and undergraduate physics textbooks, Am. J. Phys. 59 (2) (1991) 118–124. R. Scott, In the Wake of Tacoma: Suspension Bridges and the Quest for 676 Aerodynamic Stability, ASCE Press, Reston, Virginia, 2001. S.H. Strogatz, D.M. Abrams, A. McRobie, B. Eckhardt, E. Ott, Theoretical mechanics: crowd synchrony on the millennium bridge, Nature 438 (7064) (2005) 43–44. P. Dallard, A.J. Fitzpatrick, A. Flint, S. LeBourva, A. Low, R.M.R. Smith, M. Willford, The London millennium footbridge, Struct. Eng. 79 (22) (2001) 17–33. B.N. Datta, Numerical Methods for Linear Control Systems: Design and Analysis, Elsevier Academic Press, San Diego, 2004. B. Porter, R. Crossley, Modal Control: Theory and Applications, Barnes Noble, New York, 1972. B.N. Datta, Y. Saad, Arnoldi methods for large sylvester-like observer matrix equations, and an associated algorithm for partial spectrum assignment, Linear Algebra Appl. 154-156 (1991) 225–244, http://dx.doi.org/10.1016/0024-3795(91)90378-A. B.N. Datta, S. Elhay, Y. Ram, Orthogonality and partial pole assignment for the symmetric definite quadratic pencil, Linear Algebra Appl. 257 (1997) 29–48. Y. Saad, Projection and deflation methods for partial pole assignment in linear state feedback, IEEE Transactions on Automatic Control 33 (3) (1988) 290–297. B. Datta, D. Sarkissian, Partial eigenvalue assignment in linear systems: existence, uniqueness and numerical solution, In: Proceedings of the Mathematical Theory of Networks and Systems, Notre Dame, 2002. B.N. Datta, D.R. Sarkissian, Feedback control in distributed parameter gyroscopic systems: A solution of the partial eigenvalue assignment, Mech. Syst. Signal Process. 16 (1) (2002) 3–17. T.T. Soong, Active Structural Control: Theory and Practice, Longman & Scientific Technical, Essex, England, 1990. T.T. Soong, M.C. Costantinou (Eds.), Passive and Active Structural Vibration Control in Civil Engineering, Springer-Verlag, Wien GmbH, 1994. G. Housner, L. Bergman, T. Caughey, A. Chassiakos, R. Claus, S. Masri, R. Skelton, T. Soong, B. Spencer, J. Yao, Structural control: Past, present, and future, Journal of Engineering Mechanics 123 (9) (1997) 897–971, http://dx.doi.org/10.1061/(ASCE)0733-9399(1997)123:9(897). B. Spencer Jr., S. Nagarajaiah, State of the art of structural control, J. Struct. Eng. 129 (7) (2003) 845–856, http://dx.doi.org/10.1061/(ASCE)0733-9445 (2003)129:7(845). A. Preumont, Vibration Control of Active Structures: An Introduction, Springer, New York, 2011. S. Korkmaz, A review of active structural control: challenges for engineering informatics, Computers and Structures 89 (23-24) (2011) 2113–2132, http://dx.doi.org/10.1016/j.compstruc.2011.07.010. E.A. Bossanyi, Wind turbine control for load reduction, Wind Energy 6 (3) (2003) 229–244, http://dx.doi.org/10.1002/we.95. S.J. Johnson, J.P. Baker, C.P. van Dam, D. Berg, An overview of active load control techniques for wind turbines with an emphasis on microtabs, Wind Energy 13 (23) (2010) 239–253, http://dx.doi.org/10.1002/we.356. L.Y. Pao, K. Johnson, Control of wind turbines, Control Systems, IEEE 31 (2) (2011) 44–62, http://dx.doi.org/10.1109/MCS.2010.939962. A. Staino, B. Basu, S. Nielsen, Actuator control of edgewise vibrations in wind turbine blades, J. Sound Vib. 331 (6) (2012) 1233–1256, http://dx.doi.org/ 10.1016/j.jsv.2011.11.003. A. Staino, B. Basu, Dynamics and control of vibrations in wind turbines with variable rotor speed, Engineering Structures 56 (0) (2013) 58–67, http: //dx.doi.org/10.1016/j.engstruct.2013.03.014. C. Spruce, J. Turner, Tower vibration control of active stall wind turbines, IEEE Transactions on Control Systems Technology 21 (4) (2013) 1049–1066, http://dx.doi.org/10.1109/TCST.2013.2261298. D. Abramovitch, G. Franklin, A brief history of disk drive control, IEEE Control Systems 22 (3) (2002) 28–42, http://dx.doi.org/10.1109/ MCS.2002.1003997. B.M. Chen, T.H. Lee, K. Peng, V. Venkataramanan, Hard Disk Drive Servo Systems, Advances in Industrial Control, 2nd Edition, Springer, London, 2006. C. Du, L. Xie, Modeling and Control of Vibration in Mechanical Systems, CRC Press, Taylor and Francis Group, Boca Raton, FL, 2010. G.L.G. Sleijpen, H.A.V.-d. Vorst, A Jacobi–Davidson iteration method for linear eigenvalue problems, SIAM Review 42 (2) (2000) 267–293. B.N. Datta, Numerical linear algebra and applications, Society for Industrial and Applied Mathematics, Philadelphia (2010) 387–483.

S. Das et al. / Mechanical Systems and Signal Processing 72-73 (2016) 359–375

[29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71]

375

D.J. Ewins, Modal Testing: Theory, Practice and Application, John Wiley and Sons, New York, 2009. N.M.M. Maia, J.M.M. Silva, Theoretical and Experimental Modal Analysis, Research Studies Press Ltd., England, 1997. K.G. McConnell, P.S. Varoto, Vibration Testing: Theory and Practice, John Wiley and Sons, New York, 2008. S. Brahma, B. Datta, An optimization approach for minimum norm and robust partial quadratic eigenvalue assignment problems for vibrating structures, J. Sound Vib. 324 (3) (2009) 471–489. Z.-J. Bai, B.N. Datta, J. Wang, Robust and minimum norm partial quadratic eigenvalue assignment in vibrating systems: a new optimization approach, Mech. Syst. Signal Process. 24 (3) (2010) 766–783. Z.-J. Bai, M.-X. Chen, B.N. Datta, Minimum norm partial quadratic eigenvalue assignment with time delay in vibrating structures using the receptance and the system matrices, J. Sound Vib. 332 (4) (2013) 780–794. R. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer, New York, 1991. C. Soize, A nonparametric model of random uncertainties for reduced matrix models in structural dynamics, Probab. Eng. Mech. 15 (3) (2000) 277–294. V. Dragan, T. Morozan, A.-M. Stoica, Mathematical Methods in Robust Control of Linear Stochastic Systems, Springer, New York, 2006. J. Suhardjo, B. Spencer Jr., A. Kareem, Frequency domain optimal control of wind excited buildings, J. Eng. Mech. 118 (12) (1992) 2463–2481, http://dx. doi.org/10.1061/(ASCE)0733-9399(1992)118:12(2463). B. Spencer Jr., J. Suhardjo, M. Sain, Frequency domain optimal control strategies for aseismic protection, J. Eng. Mech. 120 (1) (1994) 135–158, http: //dx.doi.org/10.1061/(ASCE)0733-9399(1994)120:1(135). I.E. Kose, W. Schmitendorf, F. Jabbari, J. Yang, h1 active seismic response control using static output feedback, J. Eng. Mech. 122 (7) (1996) 651–659, http://dx.doi.org/10.1061/(ASCE)0733-9399(1996)122:7(651). L.-T. Lu, W.-L. Chiang, J.-P. Tang, M.-Y. Liu, C.-W. Chen, Active control for a benchmark building under wind excitations, J. Wind Eng. Ind. Aerodyn. 91 (4) (2003) 469–493, http://dx.doi.org/10.1016/S0167-6105(02)00431-2. B. Connor, S. Iyer, W. Leithead, M. Grimble, Control of a horizontal axis wind turbine using h infinity control, in: First IEEE Conference on Control Applications, vol. 1, 1992, pp. 117–122. http://dx.doi.org/10.1109/CCA.1992.269889 R.V. Grandhi, I. Haq, N.S. Khot, Enhanced robustness in integrated structural/control systems design, AIAA J. 29 (7) (1991) 1168–1173. K. Zhou, J.C. Doyle, K. Glover, Robust and Optimal Control, Prentice Hall, New Jersey, 1996. S.-G. Wang, H.Y. Yeh, P.N. Roschke, Robust control for structural systems with parametric and unstructured uncertainties, J. Vib. Control 7 (5) (2001) 753–772. S.-G. Wang, H.Y. Yeh, P.N. Roschke, Robust control for structural systems with structured uncertainties, Nonlinear Dyn. Syst. Theory 4 (2) (2004) 195–216. J.-C. Wu, H.-H. Chih, C.-H. Chen, A robust control method for seismic protection of civil frame building, J. Sound Vib. 294 (1–2) (2006) 314–328, http: //dx.doi.org/10.1016/j.jsv.2005.11.019. H. Du, N. Zhang, F. Naghdy, Actuator saturation control of uncertain structures with input time delay, J. Sound Vib. 330 (18–19) (2011) 4399–4412, http://dx.doi.org/10.1016/j.jsv.2011.04.025. H. Yazici, R. Guclu, I.B. Kucukdemiral, M.N. Alpaslan Parlakci, Robust delay-dependent h1 control for uncertain structural systems with actuator delay, J. Dyn. Syst. Meas. Control 134 (3) (2012) 031013-1–031013-15. A. Papoulis, S.U. Pillai, Probability, Random Variables and Stochastic Processes, Tata McGraw-Hill, New Delhi, 2002. J. Jacod, P. Protter, Probability Essentials, Springer, New York, 2004. C. Soize, Maximum entropy approach for modeling random uncertainties in transient elastodynamics, J. Acoust. Soc. Am. 109 (5) (2001) 1979–1996. S. Das, R. Ghanem, A bounded random matrix approach for stochastic upscaling, SIAM Multiscale Model. Simul. 8 (1) (2009) 296–325. S. Das, New positive-definite random matrix ensembles for stochastic dynamical system with high level of heterogeneous uncertainties, Presented at SIAM Conference on Uncertainty Quantification, Raleigh, NC, 2012. R. Ghanem, S. Das, Hybrid representations of coupled nonparametric and parametric models for dynamic systems, AIAA J. 47 (4) (2009) 1035–1044. M. Padulo, M.D. Guenov, Worst-case robust design optimization under distributional assumptions, Int. J. Numer. Methods Eng. 88 (8) (2011) 797–816, http://dx.doi.org/10.1002/nme.3203. G. Fishman, Monte Carlo: Concepts, Algorithms, and Applications, Springer, New York, 2003. C.P. Robert, G. Casella, Monte Carlo Statistical Methods, Springer, New York, 2004, 90–95. R.E. Melchers, Structural Reliability Analysis and Prediction, John Wiley and Sons, New York, 1999. T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer, New York, 2009, 520–528. MATLAB, version 8.3.0.532 (R2014a), The MathWorks Inc., Natick, Massachusetts, 2014. G. Taguchi, Introduction to Quality Engineering: Designing Quality into Products and Processes, Technical Report, Asian Productivity Organization, Tokyo, 1986. M.S. Phadke, Quality Engineering using Robust Design, Prentice-Hall, Englewood Cliffs, NJ, 1989. G. Box, S. Jones, Designing products that are robust to the environment, Total Qual. Manag. Bus. Excell. 3 (3) (1992) 265–282. J. Su, J.E. Renaud, Automatic differentiation in robust optimization, AIAA J. 35 (6) (1997) 1072–1079. G.J. Park, T.H. Lee, K.H. Lee, K.H. Hwang, Robust design: an overview, AIAA J. 44 (1) (2006) 181–191. D. Müllner, Modern hierarchical, agglomerative clustering algorithms, 2011, arXiv:1109.2378[stat.ML]. G. Carlsson, F. Mémoli, Characterization, stability and convergence of hierarchical clustering methods, J. Mach. Learn. Res. 11 (2010) 1425–1470. K. Deb, Multi-Objective Optimization using Evolutionary Algorithms, John Wiley and Sons, West Sussex, England, 2001. R.P. Santa-Ana, E. Miranda, Strength reduction factors for multi-degree-of-freedom systems, in: Proceedings of the 12-th World Conference on Earthquake Engineering, no. 1446, Auckland, New Zealand, 2000. S. Ray Chaudhuri, R. Villaverde, Effect of building nonlinearity on seismic response of nonstructural components: a parametric study, J. Struct. Eng. 134 (4) (2008) 661–670.