An Efficient Stochastic Linearisation Procedure for the Seismic Optimal Design of Passive Control Devices

The seismic response of structures is often enhanced by introducing passive control devices that can operate through the dissipation of the input energy or by modifying the dynamic characteristics of the main structure. The inherent non-linearities in the constitutive laws of some of them lead to computation difficulties and have limited the large-scale use and design of these devices. In this study, a procedure for the optimal design of multi passive control devices is proposed. The general case of linear Multi-Degree-Of-Freedom (MDOF) not-classically-damped structural systems controlled by Fluid Viscous Dampers (FVD) are investigated in a stochastic framework. The procedure consists of evaluation of the device optimal pattern by minimizing an objective function related to the dampers cost and subjected to a constraint on the structural behavior. For each configuration, the complete probabilistic characterization of the response is achieved by employing random vibration theory, Stochastic Linearisation (SL) techniques and a novel analytic model which provides closed-form PSD functions of earthquakes accelerations coherent to response spectra suggested by seismic codes. Exploiting this model, a procedure to speed up the Stochastic Linearisation technique by avoiding any numerical integration is proposed. Applications on MDOF building structures have been carried out to validate the proposed approach in terms of accuracy and reduction of the computational effort and to obtain optimal pattern of the passive control device coherently with the provisions of seismic building codes.


INTRODUCTION
The solution of a structural design problem generally requires the evaluation of a set of parameters in order to fulfill several requirements, for instance, in terms of strength, serviceability and dynamic performance of the structural system at hand. Each requirement is usually indicated as a limit state and the design problem, in other words, can be expressed as the measurement of the violation of a given set of limit states (Melchers, 1999). Since there are infinite sets of design parameters that ensure the respect of the limit states, optimization techniques may be deployed, in order to choose those parameters that satisfy the requirements and, at the same time, minimize a properly defined cost function.
Most of the established design strategies define structural safety in a deterministic way, assuming that both structural parameters and applied loads are known, even if conservative values, derived from statistical studies, are utilized. This approach leads to consider only the mean response and to neglect its dispersion. Since the uncertainties affect both structural properties and load characteristics, a more reliable analysis should be defined into a probabilistic framework. In the last decades, several methodologies have been developed for the probabilistic measure of the structural safety, in which any uncertainty about a design variable, a structural property or a load feature has to be taken into account explicitly and modeled through its respective probability density function (Melchers, 1999). Accordingly, also limit states and cost functions have to be described in terms of their probability of occurrence.
The dynamic response of structural systems subjected to seismic loads is often enhanced by introducing passive control devices, such as Fluid Viscous Dampers (FVD), Tuned Liquid Column Dampers (TLCD) or Non-linear Energy Sinks (NES). These devices operate through several different physical mechanisms by increasing the energy dissipation or by changing the dynamic characteristics of the system to which they are attached. Their use has received increasing attention in recent decades (Housner et al., 1997;Spencer and Nagarajaiah, 2003) and important applications for improving the structural performances to dynamic loads (e.g., earthquakes or winds) of new or existing building structures have been developed. The design parameters of passive control devices depend on the dynamic features of the structural system to be controlled as well as on the characteristics of the considered dynamic load. Moreover, some of the passive control devices exhibit non-linear constitutive laws and, even admitting that the primary system behaves linearly, in the controlled system a non-linear behavior may arise. This increases the computation difficulties and limits the large-scale use and design of these devices.
In fact, the probabilistic optimal design of non-linear passive control devices should be performed by means of Monte Carlo simulations, involving a very large computational effort . On the other hand, in case of linear systems, the tools of random vibration theory allow for a full probabilistic characterization of the structural response, provided that a reliable stochastic model of the seismic input is available. Aiming at this, it is common to replace the non-linear equations of motion of a system equipped with a passive device with linear equivalent ones by using well-established procedures as Stochastic Linearisation (SL) techniques (Atalik and Utku, 1976;Roberts and Spanos, 1991;Elishakoff, 2000;Alibrandi and Ricciardi, 2012). The conditions of equivalence, when passing from the non-linear to the equivalent linear system, may be established in order to preserve some response quantities in statistical sense. Since the parameters of the linear equivalent system depend implicitly on response statistics, most of SL algorithms require the solution of a system of algebraic nonlinear equations, although it has been found (Roberts and Spanos, 1991) that a simple recursive loop is adequate to simultaneously satisfy them. Hence, the computational burden required for the analyses increases, even if it is drastically reduced with respect to Monte Carlo simulations. Moreover, calculations are very often carried out only numerically, since analytical closed-form solutions for response statistics are available just for a limited class of problems (Artale et al., 2017). In Spanos and Miller (1994), formulae for the computation of the response spectral moments have been derived for the case of random excitations with band-limited white noise, Gaussian and Kanai-Tajimi seismic spectra. In order to overcome some limitations of the conventional SL, Fujimura and Der Kiureghian (2007) proposed an alternative method of linearisation, called "tail equivalent linearisation." This is a nonparametric method which consists in replacing the non-linear system with a linear one, so that the tail probability of the linear response above a specified threshold corresponds to the first order reliability approximation of the tail probability of the non-linear response above the same threshold. It provides superior accuracy for the distribution of the maximum response, especially in the tail region and it was mainly used for reliability purposes (Der Kiureghian and Fujimura, 2009).
Alternatively, Giaralis and Spanos (2010) defined a stochastic dynamic-based algorithm in order to estimate the seismic demand of bilinear hysteretic SDOF oscillators consistently with seismic code provisions that does not require non-linear numerical integrations. Their algorithm consists in the definition of a non-parametric RS-consistent stationary PSD function and in the use of the SL in order to define the so-called equivalent linear parameters (ELPs) and, finally, in the estimation of the peak inelastic response of the non-linear oscillator. Since the applicability of this latter approach was limited to relatively mild non linear response, an higher-order SL scheme was further proposed in Spanos and Giaralis (2013) for the treatment of a wider class of hysteretic constitutive laws with a resulting higher accuracy. Recently, Mitseas et al. (2018) proposed a novel stochastic dynamics framework to estimate the peak inelastic response of MDOF strongly non-linear system in a seismic context without undertaking non-linear step-by-step integrations of the response. The algorithm allows for the seismic demand estimation of MDOF systems without numerical integrations nor modal combination rules involved. Starting again from the definition of a non-parametric RS-consistent PSD function, SL is used in order to decouple the MDOF system into several SDOF oscillators characterized by ELPs. These oscillators are able to capture the peaks of the non-linear response for each DOF of interest directly in the geometric space. In order to achieve high accuracy, the algorithm requires an iterative scheme in which, for each DOF, the PSD is updated basing on the RS defined by damping modification factors.
Among the several passive control devices proposed and realized in the last decades, Fluid Viscous Dampers (FVDs) have received great attention and have been widely used in order to reduce the effects of wind or earthquakes on civil structures and to mitigate the vibrations due to shocks in mechanical equipment. Their appeal for use in seismic engineering derives from the low maintenance required, from their re-usability for subsequent earthquakes, but also from the fact that the forces exerted by the dampers and elastic forces are out of phase and the stress level in the structure is not increased by their presence. A great amount of research papers have been devoted to the optimal design and optimal placement of fluid viscous dampers to enhance the structural performance. The most relevant scientific literature in this field has been recently collected in a very comprehensive review paper  which allows for a comparison of the different strategies used to identify the optimal design configurations of FVDs. Nevertheless, it is noteworthy to recall two approaches. In the first one (Tubaldi and Kougioumtzoglou, 2015), the non-stationary stochastic response of a hysteretic SDOF system equipped with FVDs is calculated by an approximate analytical technique that makes use of a modified SL scheme and allows to consider realistic seismic excitations with time-varying frequency content. In Gidaris and Taflanidis (2015) the performance assessment and optimal design of fluid viscous dampers in a probabilistic life-cycle cost framework is discussed to obtain optimal design under different seismic scenarios.
From a mathematical perspective, the procedure presented in this work belongs to the wider class of the constrained optimization problems. In particular, the optimal pattern of the device parameters is evaluated by minimizing an objective function related to the device cost and by ensuring that a constraint function based on the structural behavior is not violated. Both cost and constraint functions are defined into a probabilistic framework. For each damper configuration, the controlled system response is evaluated by using random vibration theory in conjunction with SL technique.
In order to reduce the great computational burden involved in these problems, two novelty aspects have been here proposed and introduced. Firstly, assuming that earthquake time-histories can be modeled as samples of a stationary stochastic process, an analytical representation of earthquake load Power Spectral Density (PSD) functions consistent with seismic codes Response Spectra (RS) (Barone et al., 2015(Barone et al., , 2019 have been adopted. Secondly, an algorithm has been set up in order to obtain approximate analytical expressions of the response first spectral moments for not-classically damped Multi-degrees-of-freedom (MDOF) systems, without performing any numerical integration.
It is to remark that the more recent SL techniques for the determination of the seismic demand for hysteretic systems (Giaralis and Spanos, 2010;Spanos and Giaralis, 2013;Mitseas et al., 2018) can greatly benefit of these novelty aspects, both in terms of a significant reduction of the computational efforts in the derivation of RS-consistent PSD functions, and in accuracy, since some approximations in evaluating the response statistics can be removed without increasing the computational burden.
In the following, it will be shown how these spectral moments can be used in order to perform SL and to evaluate the optimal design of passive control systems in an efficient and accurate way. The paper is organized as in the following: in section 2 the probabilistic optimal design problem is firstly formulated in general form and then particularized for the case of n-DOF linear structural systems controlled by means of m concurrent non-linear FVDs (Di Paola et al., 2007;Di Paola and Navarra, 2009); in Section 3, in order to take advantage of the random vibration theory, the classical SL approach is described for the problem at hand and a procedure for the analytical evaluation of the response statistics, useful for SL technique, is proposed for linear n-DOF non-classically damped systems; section 4 is devoted to numerical applications in which the reduction in computational effort and the accuracy of the proposed procedure are investigated and the optimal design of FVDs is performed for a plane shear-type five-story frame and for a three-dimensional building structure. For this last application, the accuracy of the proposed procedure is further assessed in time-domain by performing a non-linear response history analysis on a set of spectrum compatible ground motion records. Finally, some conclusions are drawn. In Appendix brief details on the deriving of RS-consistent analytical PSD function and the expressions for its evaluation are provided, along with the analytical expressions that lead to the evaluation of the cross-spectral moments in modal space, once the direct spectral moments are determined.

PROBLEM FORMULATION
The equations of motion of a n-DOF linear structural system in its so-called uncontrolled state, subjected to a sample of ground motion acceleration processÜ g (t) can be expressed as: where M, C and K are the (n × n) mass, damping and stiffness matrices, respectively, U (t) is the (n × 1) vector collecting the degrees of freedom of the structural system, dots means time derivatives and τ is the (n × 1) influence vector. Aiming at reducing the dynamic response of the uncontrolled structure, it is common to introduce in the system passive control devices. The latter, operating through different physical mechanisms, generate additional damping forces or modify the dynamic characteristics of the so-called controlled system. Since many of the practically used passive control devices, such as Fluid Viscous Dampers (FVD), Non-linear Energy Sinks (NES) or Tuned Liquid Column Dampers (TLCD), present non-linear constitutive relationships, the general form of the equations of motion of the controlled system can be written as: whereÛ is the new degree-of-freedom vector, which may differ from U (t) for those categories of devices that introduce additional degrees of freedom, for example passive absorbers; g i (Ü,U,Û) and f i (t) are the total internal non-linear forces and the external applied loads, respectively, acting in the i-th degreeof-freedom direction. The dynamic behavior of the controlled system and its ability in mitigating the seismic response are governed by the design variables of the control devices, which, in turn, depend on both dynamic features of the uncontrolled system and on the characteristics of the considered dynamic load.
In general, the greater the size and cost of the used control devices, the greater reductions in response may be achieved, until allowance criteria provided by building codes are satisfied. Therefore, the problem of the optimal design of passive control devices can be properly set and, since earthquakes are random in nature, it should be tackle into a probabilistic framework, as in the following: subjected to : where φ (x, r) is the cost function to be minimized, x and r are the vectors collecting the design variables and the uncertain parameters, respectively, P s (x, r) is the survival probability,P s is the target survival probability and is the feasibility domain of the design variables. Cost function and survival probability can be expressed in terms of design variables and of the response statistics, but their analytical expressions are specific for each kind of passive control device.
Generally, the uncertain parameters r can be divided into two categories. The first one is related to the definition of the structural model (uncertainty in mass, stiffness and damping values, or model uncertainties), whereas the second one accounts for the randomness of the load. However, in case of seismic input, it is possible to neglect the uncertainties of the structural model as they lead to smaller dispersions of the response (Pinto et al., 2007). Hereinafter, the system parameters are hence considered as deterministically known and x collects only the location patterns and the design variables of the control devices.
Moreover, it is assumed that earthquake acceleration timehistories may be modeled as finite time duration samples of a zero-mean Gaussian stationary stochastic process, completely defined by the knowledge of its Power Spectral Density (PSD) function. In this case, the vector r collects the parameters of the seismic input stochastic model, obtained by using recently proposed analytical expressions of PSDs consistent with assigned RS, that cover most of the international building codes (Barone et al., 2015(Barone et al., , 2019. The Appendix reports further details of this PSD model along with the closed-form relationship between the RS parameters and the consistent PSD function ones. In recent studies (Giaralis and Spanos, 2010;Spanos and Giaralis, 2013), RS-consistent PSDs have been numerically derived and used also for stochastic analyses of non-linear yielding structures by the application of the SL. The joint use of such technique and an analytical model of the seismic PSD function could certainly improve the overall computational efficiency.
The solution of the problem expressed in Equation (3) requires, for each tentative set of the design variables, the stochastic analysis of a non-linear system, which is often performed by means of Monte Carlo simulations, involving very large computational efforts. On the other hand, in case of linear systems and a properly defined stochastic model of the seismic input, the probabilistic characterization of the structural response can be easily computed by means of the tools of stochastic analysis.
In the following, an optimal design problem is considered when FVDs are used as passive control devices. However, similar approaches may be addressed for other types of non linear passive control systems as, for instance, TLCDs (Di Matteo et al., 2014aMatteo et al., ,b, 2015 or NESs (Navarra et al., 2019b).
Viscous dampers have been widely used to mitigate the effects of wind or earthquakes on civil structures and in the shock and vibration isolation of equipment. Their appeal derives from the low maintenance required and their re-usability for subsequent earthquakes, but also from the fact that the forces exerted by the dampers are out of phase with respect to the elastic forces and do not increase the stress level in the structure. Moreover, they can be used for the protection of new constructions as well as for retrofitting purposes. Conversely, the major drawback in using viscous dampers consists in handling their non-linear force-velocity constitutive law.
The equations of motion of the controlled system in Equation (2), when m viscous dampers are concurrently deployed become: where F D (t) is the vector of the non-linear forces exerted by the viscous dampers. Although in literature there are several attempts to model the constitutive law of viscous dampers by using the theory of visco-elasticity (Schwann et al., 1988) or the fractional calculus (Makris et al., 1993), most of the manufacturers currently use a force-velocity relationship, validated by several laboratory tests, expressed as: being C d and α the characteristic parameters of the damper device,Ẏ the relative velocity at the damper ends and sign (·) the signum function. For sesmic protection purposes, values of α between 0.15 and 0.50 are used, in order to attain quite large control forces even for small relative velocities and to have almost constant output forces for large velocity values. The non-linear forces vector F D (t) can be expressed as: where Y (t) is a m-ranked vector that collects the relative displacements at the ends of each damper and all the response quantities whose statistics are useful for the computation of the cost function and of the survival probability. These quantities can be obtained as linear combinations of the degrees of freedom U (t), through the definition of a (m × n) transformation matrix R as: Since the selection of the degrees of freedom and of the response quantities of interest depends on the geometry of the problem and on the passive control devices locations, it is not possible to attain to a general expression of the matrix R. However, since FVDs are usually inserted into the frame braces, from a computational point of view, it is convenient to include in Y (t) all the inter-story drifts of the structure. The i-th element of the vector f d (Ẏ) is expressed as in Equation (5), taking into account that each damper may have, in general, different characteristic parameters C d,i and α i . Obviously, if no dampers are deployed at the i-th inter-story, a value of C d,i = 0 is set.

Cost Function and Survival Probability
In literature, several estimators of the cost function have been proposed . In Bahnasy and Lavan (2013), the use of the sum of the dampers constants C d,i at the various floors has been proposed, whereas the expected value of the sum of the peak forces of each damper has been considered in Altieri et al. (2018) and Tubaldi et al. (2016). Herein, recalling the probabilistic aspects of the proposed approach, the cost function has been defined as the sum of the characteristic values of the peak force of each damper, that represents a reliable measure of its cost, being related to its dimensioning . Hence, the cost function φ (x, r) in Equation (3) can be written as: where f dk,i is the characteristic value (at 95% percentile) of the peak force exterted by the i-th damper, σẎ i is the standard deviation of the relative velocity at the ends of the i-th damper and η k,i is the correspondent peak factor, as defined in the classical first-passage reliability theory (Vanmarcke, 1972(Vanmarcke, , 1975: In Equation (9), T s is the duration of the time window, νẎ i and qẎ i are the mean zero-crossing rate and the bandwidth factor, respectively, whose expressions can be evaluated in terms of the first spectral moments λẎ i ,j of the i-th relative velocity process: Spectral moments are defined as the geometric moments of the one-sided PSD function G X (ω) with respect to the axis ω = 0, so that the m-th order spectral moment of a generic stochastic process X(t) is given by: Spectral moments, indeed, are related to peculiar statistics of stochastic processes and allow for their characterization. For example λ X,0 and λ X,2 represent the variances of the processes X(t) andẊ(t), respectively, and other quantities as central frequency and bandwidth parameter, as well as approximate solutions for the first-passage problem can be evaluated in terms of spectral moments (Vanmarcke, 1972(Vanmarcke, , 1975(Vanmarcke, , 1976bDer Kiureghian, 1980;Di Paola and Muscolino, 1988). The survival probability P s b max,i , T s is associated to the non occurrence of crossings, into the time window T s , of the i-th maximum allowable inter-story drift b max,i , computed from the actual story height and can be evaluated as: where the risk functionα Y i b max,i is: and the quantities ν Y i , q Y i and σ 2 Y i = λ Y i ,0 are the mean zero-crossing rate, the bandwidth factor and the variance of the relative displacements, respectively. Therefore, they can be computed as in Equations (10), but with reference to the process Y i , in terms of the first spectral moments λ Y i ,j . Assuming that the failure condition is attained when only one inter-story drift exceeds the allowable value, the global survival probability P s in Equation (3) is finally obtained as:

EFFICIENT STOCHASTIC LINEARISATION TECHNIQUE
Since Equation (4) is non-linear, the evaluation of the response statistics useful for solving the optimal design problems has to be performed through Monte Carlo simulations, which require heavy computational efforts. In order to overcome this difficulty, the Stochastic Linearisation (SL) technique may be an effective tool (Atalik and Utku, 1976). The basic idea of SL is to replace the original non-linear system in Equation (4) with an equivalent linear one: in which the equivalent system matrices M (e) , C (e) and K (e) are evaluated by minimizing the difference between the two systems in statistical sense. The expressions of equivalent linear matrices for MDOF systems equipped with several kind of passive control devices, such as Tuned Liquid Column Dampers and Non-linear Energy Sinks, may be found in Navarra et al. (2019a). In the case of FVDs, it can be shown that the equivalent linear mass and stiffness matrices coincide with those of the uncontrolled system (i.e., M (e) = M and K (e) = K), while the equivalent linear damping matrix C (e) can be obtained as: in which R is the (m × n) transformation matrix, whereas the diagonal matrix D FVD can be expressed as (Di Paola et al., 2007;Di Paola and Navarra, 2009): Frontiers in Built Environment | www.frontiersin.org being Ŵ (·) the Euler Gamma function, δ ij the Kronecker delta and σẎ j = λ Y j ,2 the standard deviation of the j-th relative velocity. It is to remark that Equation (17) is derived based on a statistical equivalence in terms of damper force and by assuming the Gaussianity of the response process. Recently, other different equivalent damping formulae based on energy equivalency or by assuming other non-Gaussian probability distributions have been proposed Ricciardi, 2018, 2019). However, in all these approaches the equivalent linear system is completely defined when some response statistics are evaluated. This is, in general, achieved by using an iterative procedure. In the first iteration, it is assumed that the equivalent linear system coincides to the uncontrolled system whereas, at subsequent iterations, the estimations of equivalent linear matrices computed at the previous iteration are considered.
In conclusions, it appears that the optimal design problem into a probabilistic framework requires the evaluation in a efficient way of the first spectral moments of the response process. Furthermore, even assuming that the uncontrolled system, Equation (1), may be considered classically damped, the presence of passive control devices makes the controlled system a non-classically damped one. In this context, in order to determine the response spectral moments, two approaches may be followed. In the most common one the evaluations are carried out through numerical integrations in the geometric space, accordingly to the flowchart reported in Figure 1A. In this case, the computational burden required for the analyses increases, even if it is drastically reduced with reference to Monte Carlo simulations. In the second alternative approach, herein proposed, the computations of spectral moments, taking advantage of the analytical model of the seismic action, does not require any numerical integrations, aiming at drastically reducing the computational efforts. A flowchart of this procedure is reported in Figure 1B and these two approaches are addressed in the following subsection.

Numerical Approach for the Evaluation of Spectral Moments
The evaluation of the j-th order spectral moment matrix in terms of the response quantity of interest vector Y may be performed through numerical integration in the geometric space of the response PSD matrix: being G UU (ω) the one-sided PSD matrix of the response: in which H (ω) = [−ω 2 M (e) + iωC (e) + K (e) ] −1 is the transfer functions matrix, the asterisk denotes complex conjugate and GÜ g (ω) is the seismic load PSD function, whose expression is reported in Appendix for seismic actions consistent with several building codes RS. The spectral moments for the evaluation of the equivalent linear system, for the computation of the cost function and for the determination of the survival probability may be properly extracted from Y,j . For lightly damped and flexible systems (large period of vibration), the accuracy of Equation (18) decreases mainly due to an insufficiently long duration assumed for the stationary excitations so that the steady-state conditions are not met. In order to enhance the accuracy, a frequency-dependent corrective damping factor can be adopted (Vanmarcke, 1976a). Spanos and Giaralis (2013) provided numerical results to facilitate the selection of sufficiently long duration of stationary excitation as a function of the structural natural frequency and damping ratio.
Numerical integration of Equation (18) requires large computational efforts and may lead to inaccurate estimations of spectral moments, due to the very sharp functions involved in the case of low damping values. Since numerical integrations are required at each step of the SL, the approach above described is highly time consuming, especially when used into an optimization problem.

Analytical Evaluation of the Spectral Moments
Generally, the equivalent linear system is not classically damped and a generalized modal analysis needs to be applied. Equations of motion, as customary, can be reformulated into a set of 2n first-order differential equations: where the state-variables vector Z, the system matrix D (e) and the load vector V are defined as: The eigenvalues square roots γ i and the eigenvectors ψ i of D (e) occur in conjugate pairs and they can be collected in such a way that γ i = γ * i+n and ψ i = ψ * i+n . Moreover, due to the structure of the state-variables vector, the i-th eigenvector and the modal matrix are: where Ŵ = diag{ γ 1 , γ 2 , · · · , γ n } and = φ 1 , φ 2 , · · · , φ n are the diagonal matrix of the eigenvalues square roots and the reduced modal matrix in terms of only displacements, respectively. Lastly, the modal participation factors vector can be defined as p = −1 V. The i-th eigenvalue square root can be rewritten as γ i = −ζ i ω 0i ± i ω Di , where ω 0i , ω Di and ζ i are determined by: These parameters designate natural frequency, modal damping and damped frequency of the i-th modal oscillator (Igusa et al., 1984), respectively. Taking advantage by the analytical definition of the seismic input PSD function reported in Appendix, it is possible to perform a closed-form approximate evaluation of direct spectral moments in the modal space, without resorting to numerical integrations (Barone et al., 2019). Referring to the ith modal oscillator response process Q i (t), the expression of the j-th spectral moment λ j,Q i ,k is: when ω 0i falls within the k-th branch of the input PSD (see Equation 32), while the dimensionless quantities ϕ j,k and γ j,s are defined as follows: and the quantities e i , ω s and GÜ g (ω s ) are the parameters of the analytical model of the PSD function, whose meaning is detailed in the Appendix. Further details on the derivation and application of Equations (24) to (25) can be found in Barone et al. (2019). Di Paola and Muscolino (1988) demonstrated that the crossspectral moments of any order λ j,Q i ,Q k , if they exist for a given PSD function of the excitation, may be obtained recursively as linear combinations of the direct spectral moments. These analytical expressions are reported in the Appendix for the sake of simplicity. Furthermore, it is noted that equations from (20) to (23) account for unitary values of modal participation factors p. In this way, closed-from expressions for the determination of modal spectral moments-like those in Equation (24)-can be used straightforwardly. Conversely, when analytical expressions are not available for the problem at hand, only a limited number of numerical integrations should be performed in the modal space. Once the cross-spectral moments of the modal oscillators are obtained, the j-th order spectral moment of a set of quantities of interest in the geometric space defined by the vector Y can be computed as (Igusa et al., 1984): where the modal combination coefficients C r,ik , D r,ik , and E r,ik , which depend on the modal participating factors, can be determined as: C r,ik = a ri a rk , D r,ik = a ri c rk − a rk c ri , E r,ik = c ri c rk (27) being a ri and c ri the entries of the following matrices: whereas the matrix b = R p n , with p n = diag{p 1 , p 2 , · · · , p n }.
It is worth to note that this procedure can be generally applied for both classically and non-classically damped systems and it can be easily implemented in a computer program routine. In general, approximate evaluations of the spectral moments are achieved since analytical closed-form solutions for response statistics are available just for a limited class of problems (Spanos and Miller, 1994;Artale et al., 2017). In these cases, exact analytical evaluation of the direct spectral moments in the modal space may be used in place of the Equation (24), thus obtaining exact results. Ultimately, the evaluation of spectral moments of response quantities of interest needs only the estimation of 4n direct spectral moments in modal space, if a modal truncation is applied and onlyn ≤ n modal contributions are retained. The estimation of the spectral moments λ Y r ,j allows for the updating of the linear equivalent system, for the estimation of the cost function and for the evaluation of the survival probability, as described in the previous section 2.

NUMERICAL APPLICATIONS
In this section, aiming to show the validity of the proposed approach, the results of numerical applications are reported. Case studies of optimal design of fluid viscous dampers are addressed for a plane five-story shear-type frame and for an irregular three-dimensional building structure.
Firstly, it is to be remarked that the analytical procedure proposed in section 3.2 is based on the approximate evaluation of the direct modal spectral moments, Equations (24 and 25), by taking advantage of the definition of the probabilistic analytical model of the seismic input reported in Appendix. The same approach would also apply for other environmental hazards, i.e., wind or ocean waves, when a PSD model of the input is available. The issue of accuracy of the seismic PSD analytical model has been investigated in Barone et al. (2019), to which the reader is referred for further details. In Figure 2, however, a comparison between first order closed-form spectral moments (m = 0, 1, 2) and their numerical counterparts is shown for several values of natural periods and viscous damping ratios. Results are in very good agreement and small differences can be observed for high-order moments at large natural periods and damping ratios.
Therefore, when a specific application needs the evaluation of spectral moments in those cases in which larger errors can be anticipated, an alternative approach may consist in using numerical integration only for the 4n first direct modal spectral moments, instead of Equations (24,25), and in evaluating all other response statistics accordingly to section 3.2. In this way both computational efficiency and high accuracy may be obtained.
In all the following numerical applications. the seismic load is modeled accordingly to RS prescribed by Eurocode 8 (UNI ENV 1998:2005, 2005. In particular, it has been modeled as a process having an EC8-compatible PSD function, obtained considering the ground type A and a pseudo-stationary duration T s = 20 s. Nevertheless, similar results can be obtained for a wide range of building codes RS covered by the model adopted (Barone et al., 2015(Barone et al., , 2019. In Table 1, the parameters used to define the RS and the correspondent values in PSD analytical model, are reported. Lastly, it is worthwhile noting that, since viscous dampers are usually installed in the braces of the structural frames, the relative displacements at the damper ends constitute a subset of the interstory drifts at each story and in each direction. Therefore, it is convenient to denote with Y (t) the vector of all the inter-story drifts, to evaluate R accordingly, and to set the parameters C d,i = 0 for all the locations in which actually there are no dampers.

Accuracy and Computational Efficiency
In order to investigate the proposed procedure both in terms of accuracy and computational efficiency, a plane five-story sheartype frame, originally proposed in Takewaki (1997) and then modified in Trombetti and Silvestri (2004) has been considered. In this structural system, referred to as Model 1, both horizontal stiffness of the columns and the story mass do not vary along the height. In particular, the mass and stiffness at each floor are m i = 8 · 10 4 kg and k i = 4 · 10 7 N/m, respectively. It is also assumed a constant modal damping ratio ζ i = 0.02 and the horizontal displacements U i at each story have been chosen as degrees of freedom. The location vector is given as τ = 1 1 1 1 1 T .
In this application, the inter-story drifts are obtained as Y (t) = RU (t), where the linear transformation matrix R is: non-Gaussian SLT (EE-NG) criteria are used. Response statistics are computed for uncontrolled and controlled system and, in the latter case, for each of the aforementioned equivalence criteria. Results are reported in Table 2 in terms of computational time and of mean relative errors in the evaluation of the displacements, U i and relative velocity standard deviationsẎ i . The two mean errors are defined as follows: where the subscripts ana and num stand for analytical and numerical evaluation, respectively. Outcomes of computational times for uncontrolled system show that a great reduction of more than 90% is achieved for each evaluation of response statistics, whereas the mean errors with respect to numerical integrations, in terms of both displacements and relative velocities, are negligible.
The advantage in using analytical procedure becomes more and more evident when it is used inside a SL procedure. In fact, for the present case, eleven iterations are needed, irrespective of the equivalence criterion adopted; the computational efforts of the numerical approach increase proportionally, while no significant increments of computational time are observed for the analytical procedure. During the iterations, discrepancies due to the use of analytical procedure tend to accumulate, but the mean errors are limited below few percentage points.

Optimal Patterns and Optimal Design
In this section, the proposed analytical procedure is applied to determine the optimal pattern of fluid viscous dampers, accordingly to the minimization problem outlined in Equation (3). The minimization procedure is carried out by using the  routine fmincon in MATLAB environment that uses the interiorpoint algorithm. Following this approach, the constrained minimization is reduced to a sequence of approximate minimization problems (Byrd et al., 2000).
The probabilistic constraint in the optimization loop is defined by assuming that structure fails when the maximum inter-story drift ratio at any story and in any direction exceeds the value of 0.5%, as suggested by Eurocodes. Such a value ensures almost elastic behavior of the structure. The minimum allowable survival probabilityP s is treated as a parameter and numerical analyses are performed for different values ofP s .
With reference to the aforementioned Model 1, the optimal design of viscous dampers is performed by considering the presence of a device at each floor and assuming a constant story height h i = 3.50 m. Optimal patterns have been evaluated for several values of the dampers coefficient α, ranging from 0.15 to 1.00, and for two values of the target survival probabilityP s , 0.90 and 0.98. Figure 3 depicts the cost function of Equation (8) against α, whereas the correspondent optimal values of C d,i are shown in Figure 4.
It is worth to note that, for every value of α, the optimal pattern remains practically the same and the presence of dampers at fifth inter-story is excluded (i.e., C d,5 = 0). Moreover, the use of strongly non linear devices has proven to be beneficial as it allows for a great reduction of the maximum forces exerted by dampers, thus reducing their size and, consequently, their cost.
In Figure 5, cumulative distribution functions of the interstory drifts have been computed for probability values from 1% to 99% and compared with the maximum allowable drift (0.5% of the inter-story height). For the uncontrolled structural system (Figure 5A), i.e., when no dampers devices are used, all interstory drifts largely exceed the limit, whereas Figures 5B,C depict the configurations correspondent to the optimal deployment of viscous damper devices for α = 1.00 and α = 0.30, respectively. In the latter configurations, optimal control of FIGURE 6 | Geometrical characteristics of the six-story three-dimensional building.
structural displacements is achieved, since most of inter-story drifts reach the maximum allowable value at the same time, thus obtaining an almost uniform distribution of deformations and stresses throughout the height of the frame (Connor, 2003). The second structural model, denoted as Model 2, consists of an irregular three-dimensional six-story building structure, whose geometrical characteristics are illustrated in Figure 6.
Each floor height is 4.00 m, the first three floors have an overall dimension of 8.0 m by 16.0 m, while the top three story have plan dimension of 8.0 m by 8.0 m. All the columns are made up by HE400A profiles, whereas beams have IPE360 cross section. An uniform gravity load of q = 5.5 kN/m 2 is considered. This structure is representative of a low-tomedium rise civil steel building. Under the assumptions of rigid diaphragms at each floor and of inextensible columns, the  Frontiers in Built Environment | www.frontiersin.org dynamically significant degrees of freedom of the i-th floor can be reduced to two translations (u i and v i ) and a rotation ϕ i , for a total of 18 degrees of freedom. A constant modal damping of ζ i = 0.02, (i = 1, . . . , 6) is assumed. The structure is composed of two frames in x-direction and of three frames in y-direction. In this case, the vector Y collects all the 27 inter-story drifts (six for the frames X1, X2, Y1, and Y2 and three for the frame Y3), useful for evaluate the constraint function and the R matrix is derived accordingly.
The first columns of the Table 3 reports the modal parameters of the stiffer modal oscillators in terms of natural frequencies and of mass participating ratios for two epicentral directions of the ground motion, namely θ = 0 • (x-direction) and θ = 90 • (y-direction). For each frame, Figure 7 depicts the distribution of inter-story drifts peak values along the height computed for 25%, 50% (mean value), 75% and 90% percentiles. It can be observed that the building exhibits inter-story drift ratios that exceed the value of 0.5% for both the epicentral  Frontiers in Built Environment | www.frontiersin.org directions, even if these are larger when earthquake strikes in y-direction. Four typologies of dampers are considered for optimization purposes, all having the non-linear coefficient α = 0.3. In particular, the characteristic parameters of the dampers at the three lower floors of x-direction and y-direction frames have been denoted as C d,xl and C d,yl , respectively, whereas the correspondent parameters of the dampers at the three higher floors have been denoted as C d,xh and C d,yh , respectively.
Two analyses have been performed in order to evaluate the optimal pattern of the viscous dampers for a target survival probability ofP s = 0.90 and for θ = 0 • (OPT-X) and θ = 90 • (OPT-Y). In Table 4, the cost function and the values of the characteristic parameters of the four damper typologies are reported for θ = 0 • and θ = 90 • . It is to be noted that, irrespective of the epicentral direction, the optimal damper configurations require dampers only at lower floors, since the dynamic behavior of the building structure is dominated by the first bending mode. Table 4 also shows that the maximum interstory drift occurred at the first floor of X1 frame for θ = 0 • and at the first floor of Y1 frame for θ = 90 • . Moreover, in the last columns of Table 3 the values of equivalent linear damping at each of the first six modes are reported for uncontrolled system (UNC) and for the optimal configurations attained for θ = 0 • (OPT-X) and θ = 90 • (OPT-Y). The corrective damping ratios ζ i /ζ (Spanos and Giaralis, 2013), also reported in Table 3, show that, for the problem at hand, the stationarity assumption has negligible effects in the determination of the response statistics. Since optimal design leads to inter-story drifts close to the limit value, larger dampers are required in y-directed frames. This is obviously reflected also in equivalent linear damping values. Figure 8 depicts the distribution along height of the peak values of the inter-story drifts resulting from the optimal design for the epicentral directions θ = 0 • and θ = 90 • .
Since the proposed approach contains several approximations (mainly related to the modeling of the seismic action and on the use of SL), its effectiveness in determining the optimal configuration of dampers is further assessed in the time domain. Aiming at this, a set of seven natural recorded ground motion time-histories has been selected in order to match the EC8 RS used in the previous analyses (Iervolino et al., 2010). Individual ground motion records have been scaled and the spectrum compatibility has been achieved by imposing a ±10% tolerance with respect to the target RS in the period range between 0.10 and 2.5 s. In Table 5, all the selected records are listed, together with the main characteristics of the considered earthquakes. The results of the spectrum-compatibility check are shown in Figure 9 in terms of individual and mean RS together to the target RS and the tolerance zone.
Equations (1) and (4) have been solved by means of a stepby-step fourth-order Runge-Kutta integration scheme in order to evaluate the dynamic response of the Model 2 building subjected to the selected ground motions for the two epicentral directions θ = 0 • and θ = 90 • . Peak values of the inter-story drift have been computed for each record and their mean value has been reported in blue dashed line for comparison purposes in Figure 7 and Figure 8 for uncontrolled and optimal controlled configurations, respectively. The comparison in the uncontrolled system shows that, although stochastic analysis tends to overestimate the interstory drifts, both the methodologies are able to capture the distribution of the inter-story drift along the height of the building. Also for the optimal damper configuration cases, the time-domain analysis shows a good agreement with the results of the stochastic analysis, thus demonstrating the effectiveness of the proposed approach in the optimal design of a passive control system. However, it is to be remarked that the timedomain analysis allows to evaluate only the mean values of the response peaks, whereas the stochastic analysis provides the full probabilistic characterization of the response.

CONCLUSIONS
In this paper, a methodology for the optimal design of passive control systems into a probabilistic framework has been described. In this way, uncertainties and response dispersion have been taken into account and, at the same time, a seismic analysis consistent with the Eurocode 8 response spectra has been carried out. The following conclusions can be drawn: • A general form of Stochastic Linearisation approach for the cases of MDOF structural systems controlled by multi concurrent passive control devices has been presented and an iterative efficient solution, able to avoid any use of numerical integrations, has been proposed for the case of seismic loads. • Herein, an optimal design procedure for fluid viscous dampers has been formulated. For this purpose, a cost function has been used as objective function to minimize. This has been defined as the sum of the characteristic values of the forces exerted by the dampers themselves. • Taking advantage of an analytical model for the computation of response spectra compatible PSD functions as well as of complex modal analysis in state-space, analytical evaluations of the first spectral moments of a set of response quantities of interest have been derived for not-classically damped systems. • The design problem has been posed in terms of survival probability considering that the failure of the structural system is identified with the over-crossing of a maximum allowable relative displacement. • The validity of the proposed approach has been investigated in terms of accuracy and computational efficiency and two applications have been presented. The first one is related to a plane shear-type five-story frame, while the second one deals with a six-story three-dimensional building structure. Optimal damper configurations and cumulative distribution functions of inter-story drifts have been evaluated and it is proved that the proposed procedure leads to an accurate evaluation of response statistics. • The proposed procedure allows for a dramatic reduction of the computational time, especially when is used for those problems that make intensive use of stochastic linearisation technique, as optimal design of passive control device. • The effectiveness of the proposed approach has been assessed against the results of time-domain analyses by using a set of seven scaled spectrum-compatible natural ground motion records. • The most recent SL techniques for the hysteretic systems seismic demand evaluation can benefit of the two novel aspects herein presented, namely the analytical model of PSD functions and the algorithm for evaluating spectral moments, in terms of both computational efficiency and accuracy.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
All the authors contributed in an equal manner to the study of the relevant literature, to the derivation of the proposed approach, to the preparation of software routines for the numerical applications and to writing and editing of the manuscript.

Stochastic model of the seismic action consistent with Response Spectra
The procedure for the optimal design of passive control devices into a probabilistic framework requires the modeling of the seismic ground motion as a stochastic process. Generally, international building codes define the seismic action by means of pseudo-acceleration elastic uniform hazard spectra associated with the peak response of linear single-degree-of-freedom systems having viscous damping. A general class of RS can be expressed as in the following four-branches expression: where T is the natural period of the SDOF system, S 0 is the peak ground acceleration, a is the dynamic amplification factor, T 1 , T 2 and T 3 are the periods that define the various branches, and k 1 and k 2 are shape factors. In this paper we refer only to RS with k 1 = 1 and k 2 = 2 and in Figure A1A a qualitative representation of RS is depicted. Several building codes allow for an alternative representation of the seismic ground motion by means of artificial accelerograms of nominal duration T s that can be generated as samples of a zero-mean Gaussian stationary process, fully characterized by its one-sided PSD function GÜ g (ω).
In this work, the RS-consistent PSD was intended as an alternative conventional way to define the seismic action or, in other word, as a mathematical tool to conveniently represent the seismic action and to be used into the linear stochastic dynamic context. With this in mind, the RS-consistent PSD function is determined by solving an inverse stochastic dynamic problem in order to produces the same effects of the target RS, hence in Equation (A2) a value of ζ = 0.05 has been used. Once the PSD function was determined, it completely defines the input process from a probabilistic point of view, irrespective of the characteristics of the superimposed structure (linearity, damping values and so on). It is remarked that in some recent SL approaches (Mitseas et al., 2018) the PSD is updated at each iteration to be consistent with a RS in which damping modification factors are applied.
Seismic codes do not define the process, but require, instead, that it has to be compatible with an assigned RS, by providing the compatibility conditions. For instance, following the provision of Eurocode 8 (UNI ENV 1998:2005, 2005, a ground acceleration PSD function GÜ g (ω) is considered compatible with an assigned acceleration RS, S a (T), if a SDOF system with an assigned damping ratio (usually ζ 0 = 0.05), subjected to accelerogram samples generated from GÜ g (ω), experiences into a time window of the nominal duration T s of the earthquake an average absolute peak acceleration larger than 90% of S a (T) for each value of the natural period T. If the ground motion PSD was known, the corresponding RS could be easily obtained by stochastic analysis. However, the inverse problem (i.e. determining the PSD function corresponding to an assigned RS) is not easy to solve. An approximate recursive solution for this problem has been provided in (Cacciola et al., 2004) in order to obtain an estimate of the PSD function GÜ g (ω) compatible with the assigned RS: where the parameter γ = 4ζ / (π − 4ζ ) and the peak factor η U (ω, ζ ) is computed for 50% probability of nonexceedance. In this choice, a further little approximation is introduced, i.e. by confusing the mean value with the median.
Recently, an analytical model of PSD function, compatible with RS building code defined by Equation (A1), has been proposed (Barone et al., 2015(Barone et al., , 2019 as: where ω i = 2π/T i (i = 1, 2, 3) is the circular frequencies corresponding to the periods T i (i = 1, 2, 3) and G 0 represents the peak value of the PSD function that occurs at the frequency ω = ω 2 . The proposed model, whose graphical representation is reported in Figure A1B, depends on only five parameters, namely G 0 and the four exponents e 1 , ..., e 4 .