Tradeoff analysis of non-linear dynamical systems under stochastic excitation

Probabilistic Engineering Mechanics 24 (2009) 585–599

Contents lists available at ScienceDirect

Probabilistic Engineering Mechanics journal homepage: www.elsevier.com/locate/probengmech

Tradeoff analysis of non-linear dynamical systems under stochastic excitation H.A. Jensen Department of Civil Engineering, Santa Maria University, Casilla 110V, Valparaiso, Chile

article

info

Article history: Received 2 September 2008 Received in revised form 20 March 2009 Accepted 6 April 2009 Available online 15 April 2009 Keywords: Global approximations Multi-objective optimization Pareto points Sensitivity analysis Statistical linearization Stochastic process Tradeoff analysis

abstract An efficient methodology to carry out multi-objective optimization of non-linear structural systems under stochastic excitation is presented. Specifically, an efficient determination of particular Pareto or noninferior solutions is implemented. Pareto solutions are obtained by compromise programming which is based on the minimization of the distance between the point that contains the individual optima of each of the objective functions and the Pareto set. The response of the structural system is characterized in terms of the first two statistical moments of the response process, i.e. the mean and variance. An efficient sensitivity analysis of non-inferior solutions with respect to the design variables becomes possible with the proposed formulation. Such information is used for decision making and tradeoff analysis. The compromise programming problem is solved by an efficient procedure that combines a local statistical linearization approach, modal analysis, global approximation concepts, and a sequential optimization scheme. Numerical results show that the total number of stochastic analyses required during the multiobjective optimization process is in general very small. Hence, different compromise solutions including the design that best represents the outcome that the designer considers potentially satisfactory are obtained in an efficient manner. In this way, the analyst can conduct a decision-making analysis through an efficient interactive procedure. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction The design and optimization of complex engineering problems usually involve several design goals which are required to be maximized or minimized while satisfying a number of design constraints. These design goals are potentially conflicting requirements reflecting technical and economical performances of a given system design. To accommodate these conflicting design goals and to explore the design options, the optimization problem is formulated in terms of a multi-objective optimization procedure [1]. In fact, multi-objective optimization has proven to be an efficient tool for solving such conflicts in engineering design. An important class of solutions to the multi-objective problem is said to belong to the Pareto front (set of Pareto optimal solutions). Each solution or alternative design comprising the front is understood to be Pareto optimal which means there are no other design alternatives for which all objectives are better. Mathematically, the Pareto optimality is defined as follows [2–5]: A vector of design variables {θ ∗ } is a Pareto optimum if and only if, there is no other design {θ } such that fj ({θ }) ≤ fj ({θ ∗ }) for all j = 1, . . . , m and fj ({θ}) < fj ({θ ∗ }) for at least one j, 1 ≤ j ≤ m, where fj , j = 1, . . . , m represent the objective functions involved in the problem. In other words, for any other design different than {θ ∗ }

E-mail address: [email protected]. 0266-8920/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.probengmech.2009.04.003

at least one objective function increases. The selection of one of the Pareto design alternative over another is a matter of decisionmaker preference. The Pareto front includes the individual optima of each objective function and it is in general a non-linear surface that may contain discontinuities. Classical methods for generating the Pareto optimal set combine the objectives into a single parameterized function. These methods involve assigning a weight to each objective function. Several optimization runs with different settings are performed in order to generate the Pareto optimal set. The results of the identification of a particular design depend on the set of weights used in the analysis. The single objective problem is computationally attractive since conventional minimization algorithms can be applied to solve the problem [6–8]. The main limitation of weighted sum methods is that they do not yield solutions that lie in non-convex regions of the feasible design space. In practice this means that weighted sum methods could miss potentially preferable designs. Additional procedures available for the generation of compromise solutions with various capabilities and limitations, include goal programming, physical programming, and interactive methods [9,10]. Another methods that have also been used for solving multiobjective optimization problems are based on genetic algorithms [11–14]. Genetic algorithms process a set of promising solutions simultaneously and therefore they are capable of capturing several points along the Pareto front. These algorithms are based on

586

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

an arbitrary initialized population of search points in the parameter space, which by means of selection, mutation and crossover evolves towards better regions in the search space. These methods are able to search a greater design space than other optimization procedures but they are associated with high computational costs. In general, the determination of the entire Pareto front is prohibitively expensive for complex structural systems. This is the case of non-linear dynamical systems under stochastic loading since the process may required a considerable number of stochastic analysis calls. The aim of this work is restricted to obtain in an efficient manner, compromise solutions that best represent the outcome that the designer considers potentially satisfactory and the corresponding sensitivity information of such solutions. This type of information provides useful information to the designer or analyst for decision making and tradeoff analysis, since it allows an efficient design exploration around specific Pareto points. The compromise designs are obtained by solving a sequence of explicit optimization problems whose solutions can be obtained in an efficient manner. The response of the structural system is characterized in terms of the first two statistical moments of the response process, i.e. the mean and variance. The evaluation of the second moment characteristics of the response process for a given design is carried out by an efficient procedure that combines an orthogonal series expansion of the stochastic processes involved in the problem and a local linearization technique. The efficiency of the procedure is further increased by considering analytical approximations of the objective and constraint functions in terms of the design variables during the optimization process. The analytical approximations are constructed by combining a mixed linearization approach with a stochastic response sensitivity analysis. Two example problems, including a large non-linear finite element model, are presented to illustrate the applicability of the proposed procedure.

3. Compromise programming The determination of a set of non-dominated solutions (Pareto optimum solutions or non-inferior solutions) of the above multiobjective optimization problem achieves a compromise among different objective functions. As previously pointed out there are a number of procedures available for the generation of compromise solutions with various capabilities and limitations. The objective of this work is to obtain in an efficient manner particular compromise solutions. To this end, the set of compromise solutions is defined by means of the so-called ideal point, F (id) , which contains the individual optima of each of the objective functions. The approach, which is called compromise programming (CP) is based on the minimization of the distance between the ideal point and the Pareto set [1]. In particular, the Tchebyshev norm [15] is used to measure such distance in the present implementation. The advantage of the Tchebyshev norm approach is that efficient solutions are produced even when the Pareto set is non-convex. The corresponding compromise programming problem to the multi-objective optimization problem (1) can be written as Min d + α

ωi σfi ({θ }) +

i =1

ndf X

# ω

d d i fi

({θ })

i=1

d≥

d≥

µfi ({θ}) − µ(fiid) µ(fias) − µ(fiid) σfi ({θ}) − σf(i id) σf(i as) − σf(i id) fid ({θ}) − fi fi

d(as)

d(id)

− fid(id)

,

i = 1, . . . , nsf

,

i = 1, . . . , nsf

,

i = 1, . . . , ndf

gj µri ({θ }), σri ({θ }), rld ({θ}) ≤ 0




i = 1, . . . , nsr , l = 1, . . . , ndr , j = 1, . . . , nc , {θ} ∈ Θ

(2)

where α is a small positive number, d is a positive auxiliary µ variable, and the weights ωi , ωiσ and ωid are defined as

subject to

µ

ωi =

gj µri ({θ }), σri ({θ }), rld ({θ }) ≤ 0




i = 1, . . . , nsr , l = 1, . . . , ndr , j = 1, . . . , nc , {θ} ∈ Θ

ωi µfi ({θ}) +

σ

subject to

d≥

n {F ({θ})} = µf1 ({θ }), . . . , µfnsf ({θ }), σf1 ({θ}), . . . , o σfnsf ({θ}), f1d ({θ }), . . . , fnddf ({θ })

nsf X

µ

i=1

2. Problem formulation Consider a multi-objective optimization problem defined as the identification of a vector {θ}, θi , i = 1, . . . , nd , of design variables to minimize a vector of objective functions

"n sf X

(1)

where µfi ({θ}) and σfi ({θ})are the maximum values of the mean and standard deviation of the objective function fi , respectively, fld ({θ }) is a deterministic objective function, µri ({θ }) and σri ({θ}) are the maximum values of the mean and standard deviation of the response function ri , respectively, rld ({θ}) is a deterministic response function, gj , j = 1, . . . , nc are the constraint functions, and Θ is the set that contains the side constraints for the vector of design variables. The deterministic objective functions fld ({θ }), l = 1, . . . , ndf and the deterministic response functions rld ({θ }), l = 1, . . . , ndr are assumed to be explicit functions of the design variables and they are related to design requirements such as structural weight, geometric conditions, material cost components, etc. The constraint functions gj , j = 1, . . . , nc are assumed to be explicit functions of the quantities µri ({θ}), i = 1, . . . , nsr , σri ({θ }), i = 1, . . . , nsr , and rld ({θ }), l = 1, . . . , ndr . A simple linear combination of these quantities, which is the case in most practical applications, is considered here.

ω = d i

1

µ(fias) − µ(fiid)

ωiσ =

,

1

σf(i as) − σf(i id)

1 fi

d(as)

(3)

− fid(id)

(as)

where µfi

,

> µ(fiid) , σf(i as) > σf(i id) , and fid(as) > fid(id) . The (as)

(as)

d(as)

components µfi , σfi , and fi are called the aspiration levels and they are defined by the designer or analyst prior to solving the compromise programming problem. These aspiration levels represent the outcomes that the designer feels adequate, and they can only be achieved when they are on the compromise set. In other cases, a Pareto point that best represents the aspiration levels is obtained. Fig. 1 illustrates this concept for the case of two objective functions, i.e. F1 and F2 . It is seen that if the vector of aspiration levels {F (as) } is not on the Pareto front, compromise programming gives the Pareto point {F ∗ } that best represent such aspiration levels. This point is obtained by minimizing the distance between the ideal point and the Pareto curve in the direction defined by the ideal point and the vector of aspiration levels.

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

F1

587

FEASIBLE SPACE {F(as)} Aspiration Level

{F(*)} Pareto Point PARETO CURVE

{F(id)} Ideal Point

F2 Fig. 1. Pareto point {F (∗) } that best represents the aspiration level {F (as) }.

4. Pareto sensitivity analysis The sensitivity analysis at a given Pareto solution provides the variation in one objective given the variation in another objective on the Pareto surface in a given direction. In particular, the sensitivity along the feasible descent direction of each of the objectives in the objective function space is considered here. This information can be used for updating the aspiration levels in the context of compromise programming, generating first or second-order Pareto surface approximations, conducting a tradeoff analysis and decision making, and checking the robustness of the current design. Let {θ ∗ } be the design corresponding to the noninferior solution that best represents the designer aspiration level, and J the active constraint set defined as J = {j : gj (µri ({θ ∗ }), σri ({θ ∗ }), rld ({θ ∗ })) = 0, j = 1, . . . , nc },

J = {j : − ≤ gj (µri ({θ }), σri ({θ }), ∗

rld

({θ })) ≤ 0, j = 1, . . . , nc } (5) ∗

in which  is a user-defined scalar which allows constraints near the limit to be included in the active set. If Fi represents one of the objective functions of the stochastic multi-objective optimization problem, i.e., µfi , σfi , i = 1, . . . , nsf , or fid , i = 1, . . . , ndf , then the feasible direction {dFi } with the greatest improvement of the objective Fi is obtained by projecting the gradient vector −∇θ Fi ({θ ∗ }) onto the projection matrix [P ] of the active constraint set, that is,

{dFi } = −[P ] ∇θ Fi ({θ ∗ })

(6)

where the projection matrix is given by

[P ] = [I ] − [B]T ([B][B]T )−1 [B]

(7)

and where the rows of the matrix [B] are the gradients of the active constraints (constraints in J and the active bounds) [16]. For example, if there are no active bounds the matrix [B] is defined as

 .  .    [B] = ∇ gjT  ,  .  . 

j ∈ J.

is dependency, it is necessary to modify the mathematical model by deleting the redundant constraints [17]. The Pareto sensitivity measure δ Fi /δ Fj , i.e., the variation in one objective given the variation in another objective, along the feasible descent direction of the objective Fj is given by the projection of {dFi } onto the direction {dFj }, i.e.

δ Fi /δ Fj =

(4)

where in actual computation the active set is implemented as ∗

Fig. 2. Direction {dFi } with the greatest improvement of the objective Fi in the direction of the tangent plane defined by the active constraints.

(8)

The vector {dFi } corresponds to the projection of the steepest descent direction of the objective Fi onto the tangent plane defined by the active constraints at {θ ∗ } (see Fig. 2 for a schematic illustration of the vector {dFi }). In the previous derivation it has been assumed that the current Pareto point {θ ∗ } is a regular point, that is, the rows of the matrix [B] are linearly independent. If there

{dFi }T {dFj } . {dFj }T {dFj }

(9)

This equation represents the tradeoff information of the objective Fi along the feasible descent direction of the objective Fj . In a similar manner a Pareto sensitivity analysis along the feasible descent direction of each of the objective functions can be conducted. This sensitivity information can be represented in matrix form, called the tradeoff matrix, as



1

 δ F /δ F [T ] =  1 2 ... δ F1 /δ Fm

δ F2 /δ F1 1

... δ F2 /δ Fm

... ... ... ...

 δ Fm /δ F1 δ Fm /δ F2  ... 

(10)

1

where m is the number of objective functions, in this case, m = 2nsf + ndf (see Eq. (1)). It is noted that the sensitivities predicted from the previous analysis are only valid if the active constraints remain unchanged for a small move in the feasible descent directions of the objective functions along the Pareto surface from the current non-inferior solution. The sensitivity information given by the tradeoff matrix requires the sensitivity evaluation of the quantities µfi ({θ }) and σfi ({θ }), i = 1, . . . , nsf , and fid ({θ}), i = 1, . . . , ndf (gradients of the objective functions), and µri ({θ }) and σri ({θ }), i = 1, . . . , nsr , and rid ({θ }) i = 1, . . . , ndr (definition of the matrix [B]). The evaluation of such sensitivities is discussed in subsequent sections. From the tradeoff matrix is clear that the tradeoff between any two objective functions exists only if the corresponding off-diagonal term is negative, and the magnitude of the tradeoff is given by its absolute value. On the other hand, a positive off-diagonal term indicates that there is no tradeoff between the corresponding objectives, i.e, both of them can be improved simultaneously. Note that at least one of the elements in any given row has to be negative for a given design to be a local Pareto optimum (necessary condition).

588

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

5. Local pareto surface approximation The sensitivity information previously defined can also be used to find the normal to the Pareto surface, which defines the tangent hyperplane to the Pareto surface at the current compromise design. This representation is useful in decision making because it provides an inclusive local information about the Pareto surface and facilitates the designer to explore other compromise designs. At this point it is necessary to mention the limitations and conditions under which a local first-order approximation is valid. In general the Pareto surface can be highly non-linear and may have discontinuities. Then, for smooth approximations to work the active constraint set changes should be mild, and for them to be valid in the entire vicinity of the current point there should be no nearby jumps (first-order discontinuities) in the Pareto surface. The mathematical description of the first-order approximation of the Pareto surface can be written as

X δ Fm

Fm = Fm +

i=1

δ Fi

7. Numerical implementation 7.1. Approximate stochastic responses

m−1

∗

new set of aspiration levels should be defined. In realistic problems the number of objectives can be large and the designer might not be able to set new aspiration values for each of the objective functions, but only for some of them. In this regard, the sensitivity analysis can be used to generate aspiration levels for some objectives based on the ones set for the other objectives (see Eq. (11)). These new aspiration levels can be used as approximate new Pareto solutions, or they can be employed to generate a new non-inferior solution by solving the corresponding compromise optimization problem. In summary, the sensitivity information presented in previous sections provides the designer with a formal means for efficient exploration around a given Pareto design and for decision making and tradeoff analysis.

(Fi − Fi ) ∗

(11)

where Fi∗ = Fi ({θ ∗ }), i = 1, . . . , m is the current Pareto point in the objective function space, δ Fm /δ Fi is the change in the objective Fm for a given change in Fi while the other objectives are kept constant. Eq. (11) represents one of the objective functions in terms of the others m − 1 objectives. The procedure for determining the sensitivity measure δ Fm /δ Fi in Eq. (11) (only the objectives Fm and Fi change) is similar as the one presented before except that the matrix [B] is modified to include not only the gradients of the active constraints but also the gradients of each of the objective functions that are kept constant. A second-order approximation to the Pareto surface can be computed in a similar manner. In this case, additional Pareto points can be generated by using a projection method that uses the gradient information available at the current Pareto design [10]. On the other hand, Pareto data around the current design point can also be obtained by generating Pareto points from aspiration levels in the neighborhood of the present level chosen by the designer. This data can then be used to approximate the Pareto surface (first or second-order approximations) in the vicinity of the current Pareto point.

The computational cost problem involved in the solution of the compromise optimization problem (2) is addressed by the use of approximate stochastic analyses during portion of the optimization process [18]. A stochastic analysis block is first used to analyze an initial design and then to generate information that allows the construction of objective and constraint approximations. Following an optimization based on an approximate analysis, an exact stochastic analysis is performed in the design point obtained by the approximate optimization so that a new approximation for the objective function and constraints can be constructed. The process is repeated until convergence is achieved, typically measured by the magnitude of changes in the objective function. In the proposed formulation, the first and second-order statistics of the stochastic performance functions (objectives fi , i = 1, . . . , nsf , and response functions ri , i = 1, . . . , nsr ) are represented by using approximate functions. Let τ ({θ }) be the maximum value of the first or second-order statistic of a generic performance function, and {θ 0 }, θi0 , i = 1, . . . , nd , a point in the design space. The function τ ({θ}) is approximated globally about the point {θ 0 } as

τ¯ ({θ }) = τ ({θ 0 }) +

nd X ∂τ ({θ 0 }) i=1

6. Tradeoff analysis If the current design is Pareto optimum, i.e., obtained by compromise programming, there is no other feasible design that improves all the objective functions. Therefore, if the designer wants to improve some of the objective functions it can only be done at the expenses of other objective functions. The tradeoff analysis addresses the issue of how much one has to relax some objectives for compensating the improvement of other objectives. In complex design problems it is often difficult to know how much to change a given objective to obtain a satisfactory design. The approximate analysis presented earlier provides Pareto surface exploration capabilities which can be used to guide the designer in setting the amount of relaxation or improvement of a given objective function in the vicinity of the current non-inferior solution. The Pareto sensitivity analysis and the local first or secondorder Pareto surface approximation can also be used to check the consistency in the aspiration levels set by the designer. This is important when the aspiration levels are set without full knowledge of the Pareto surface in the vicinity of the current design. In that case, the aspiration levels may be inconsistent or unattainable. It is emphasized that the only aspiration levels that are attainable are the ones that are on the Pareto surface. On the other hand, if the designer is not satisfied with the current design, a

∂θi

B(θi , θi0 ),

(12)

where B(θi , θi0 ) is an operator that defines the type of approximation for the design variable θi . In particular, the following operator is implemented

B(θi , θi0 ) =

   θi − θi0 θ0    i (θi − θi0 ), θi

∂τ ({θ 0 }) ≥0 ∂θi 0 ∂τ ({θ }) if < 0. ∂θi

if

(13)

This mixed linearization is usually called convex approximation since τ¯ ({θ}) is a convex function [19,20]. An attractive property of this approximation is that it yields the most conservative approximation among all possible combinations of direct/reciprocal variables. This property is important since it is the basis of the optimization scheme proposed in the sequel for solving the compromise optimization problem (2). Using the mixed linearization, the following explicit compromise optimization problem is generated: Min d + α

"n sf X

µ

ωi µ ¯ fi ({θ }) +

i=1

nsf X

σ

ωi σ¯ fi ({θ }) +

i=1

subject to µ

d ≥ ωi



 µ ¯ fi ({θ }) − µ(fiid) ,

i = 1, . . . , nsf

ndf X i=1

# ω

d d i fi

({θ })

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599



(id)



d(id)

d ≥ ωiσ σ¯ fi ({θ }) − σfi d ≥ ωid fdi ({θ}) − fi



,



,

i = 1, . . . , nsf i = 1, . . . , ndf

¯ ri ({θ}), σ¯ ri ({θ }), rld ({θ }) ≤ 0 gj µ


i = 1, . . . , nsr , l = 1, . . . , ndr , j = 1, . . . , nc , {θ } ∈ Θ

(14)

where nd X ∂µf ({θ 0 })

µ ¯ fi ({θ }) = µfi ({θ 0 }) +

i

∂θi

i=1

σ¯ fi ({θ }) = σfi ({θ 0 }) +

nd X ∂σf ({θ 0 }) i

∂θi

i =1

µ ¯ ri ({θ}) = µri ({θ 0 }) +

B(θi , θi0 )

nd X ∂µr ({θ 0 }) i

∂θi

i =1

σ¯ ri ({θ }) = σi ({θ 0 }) +

B(θi , θi0 )

nd X ∂σr ({θ 0 }) i

i =1

∂θi

B(θi , θi0 )

589

optimization problem is constructed about a feasible design, the feasible domain corresponding to the approximate problem tends to be inside the feasible domain of the original problem (2). As a result, the process has the tendency to generate a sequence of steadily improving feasible designs. It is noted that the above statement is not general in an absolute sense, since it is not known whether the approximations are more conservative than the exact functions, which are unknown. In this regard convex approximation of highly non-convex quantities may not be appropriate in terms of the desired properties of the explicit compromise optimization problems, i.e., conservative feasible domains. Finally, it is noted that even though the convergence of the proposed optimization scheme is not guaranteed from a strict mathematical point of view, the scheme has been found useful in a large number of structural optimization applications, including the cases considered in the context of this study. 8. Application to dynamical systems

B(θi , θi0 )

(15)

are the approximate quantities and all other terms have been previously defined. Note that in the above optimization problem the objective and constraint functions are explicit functions of the design variables {θ }. The sensitivities of the maximum value of the first two statistical moments of the stochastic performance functions are obtained by performing a stochastic response sensitivity analysis. Such analysis is introduced in subsequent sections. 7.2. Solution of compromise optimization problem As previously pointed out the solution of the compromise optimization problem (2) is obtained by transforming it into a sequence of sub-compromise optimization problems [21,18]. At the beginning of each design cycle the system is analyzed to construct an explicit approximate optimization problem of the form given in Eq. (14). The approximations are constructed about the current design point {θ k }, θik , k = 0, 1, . . . , i = 1, . . . , nd , where k indicates the present cycle in the overall optimization process. The characterization of the approximate optimization problem requires a stochastic response sensitivity analysis at the current design point. The explicit approximate optimization problem is solved by any non-linear optimization algorithm to obtain an optimal solution {θ ∗ }. The solution can be obtained in an efficient manner due to the explicitness of the quantities involved in the problem in terms of the design variables. The new point in the design space {θ ∗ } is used as the current design point for the next optimization cycle, that is, {θ k+1 } = {θ ∗ }. The design {θ k+1 } is the point where the compromise optimization problem is approximated in the next cycle. The process is continued until some convergence criterion is satisfied. In particular, a convergence criterion defined in terms of the relative change of the auxiliary variable d and the objectives µfi , σfi and fid between two consecutive iterations is adopted in this implementation. It is noted that the characterization of the approximate stochastic responses (Eq. (12)), and hence the construction of the explicit compromise optimization problem (14) requires just one stochastic analysis (sensitivity analysis) at a given design point. Similarly, the Pareto sensitivity analysis needs only one stochastic sensitivity analysis at the given or current Pareto optimal point. Due to the convex linearization defined in (12), the quantities involved in the approximate optimization problem (14) are locally conservative. This means that they tend to overestimate the values of the true functions. Then, if the approximate compromise

8.1. Mechanical modeling The multi-objective optimization formulation presented in previous sections is applied to a class of non-linear structural systems under stochastic excitation. The displacement response {u(t )} of the structural system is characterized by the equation of motion

[M ]{¨u(t )} + [C ]{˙u(t )} + {p(t )} = {b}a(t ),

(16)

where {p(t )} denotes the vector of non-linear restoring forces, {˙u(t )} the velocity vector, {¨u(t )} the acceleration vector, [M ] the mass matrix, [C ] the damping matrix, and {b} the vector that couples the stochastic process a(t ) to the degrees of freedom of the structure. The non-linear restoring force vector {p(t )} is assumed to be a function of {u(t )}, {˙u(t )}, and a set of auxiliary variables {z (t )}. It is written as {p(t )} = [K ]{u(t )} − [Q ]{z (t )} where [K ] is the linear stiffness matrix, and [Q ] is a constant matrix. Such characterization permits to represent several non-linear models of practical relevance, including hysteretic systems. The evolution of the set of auxiliary variables satisfies a non-linear differential equation of the form

{˙z (t )} = {h({u(t )}, {˙u(t )}, {z (t )})} ,

(17)

where {h(·, ·, ·)} is a non-linear vector function. The excitation a(t ) is assumed to be a continuous Gaussian second-order stochastic process, and it is expressed by the Karhunen–Loève (K–L) expansion as [22,23] a(t ) ≈ a0 (t ) +

nkl X

aj (t )ξj

(18)

j =1

where ξj , j = 1, . . . , nkl , are independent, identically distributed standard Gaussian random variables, a0 (t ) and aj (t ) denote the mean function and the jth K–L component of a(t ), respectively, and nkl is the order of truncation of the series expansion. The above representation of the stochastic process is quite general since it permits to describe any kind of Gaussian excitation including white noise, colored excitation, non-zero mean excitation, and nonstationary excitation. In addition, this representation can be easily obtained if measurements of the excitation process are available. 8.2. Stochastic response evaluation The stochastic performance functions involved in the compromise optimization problem, i.e. fi , i = 1, . . . , nsf (objective functions), and ri , i = 1, . . . , nsr (response functions) are given in terms

590

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

of the response process of the mechanical model characterized by Eqs. (16) and (17). The second-moment characterization of these quantities is obtained by applying a statistical equivalent linearization approach [23,24]. In particular, an equivalent linearization technique based on the Gaussian assumption of the response process is used in the present implementation. It is noted that the constraints on applying this technique are similar to the constraints on all other equivalent linearization based methods. Then, depending on the level and type of non-linearity equivalent statistical linearization may not be suitable for obtaining accurate estimates of second-order statistics. Numerical validations have shown that the linearization technique used in this formulation gives results with sufficient accuracy for the type of non-linearity considered in the numerical examples [25]. The computational complexity to linearize the non-linear vector function {h(·, ·, ·)} (Eq. (17)) grows significantly with increasing number of non-linear components. The linearization of the non-linear quantities is performed for each non-linear component individually improving in this manner the efficiency of the linearization procedure [26]. Assembling the linearized equations of all non-linear elements leads to a linear firstorder differential equation for the temporal evolution of the set of variables {z (t )}. The equation can be written in terms of the structural system state, expressed by {u(t )}, {˙u(t )} and {z (t )} as

{˙z (t )} = {µ{z (t )} } + [Υ1 (t )] {u(t )} − {µ{u(t )} }
+ [Υ2 (t )] {˙u(t )} − {µ{˙u(t )} }
+ [Υ3 (t )] {z (t )} − {µ{z (t )} } {µ{z (t )} } = E ({h({u(t )}, {˙u(t )}, {z (t )})})




(19)

where the matrices [Υ1 (t )], [Υ2 (t )] and [Υ3 (t )] contain the linearization coefficients with {µ{u(t )} } = E ({u(t )}), {µ{˙u(t )} } = E ({˙u(t )}), and {µ{z (t )} } = E ({z (t )}). The linearization coefficients contained in [Υ1 (t )], [Υ2 (t )] and [Υ3 (t )], and the mean vector {µ{z (t )} } depend on the system response, which in turn depends on the linearization coefficients. Based on the superposition principle and the representation of the excitation as a sum of independent loading components according to (18), the responses {u(t )}, {˙u(t )} and {z (t )} are represented in the same manner, i.e.

{u(t )} ≈ {u(t )}0 +

The vectors of the K–L representation of the response process can be computed by solving Eqs. (16) and (19). Because Eq. (19) involves statistical equivalent linearization, it depends on the second moment characteristics of the response quantities. Therefore, the vectors {u(t )}0 , {˙u(t )}0 , {z (t )}0 , {u(t )}j , {˙u(t )}j and {z (t )}j , j = 1, . . . , nkl can not be calculated directly, but in an iterative manner. First, the representation of the response process given in (20) is substituted in the equation of motion (16) and the linearized equation (19). Next, the equation of motion expressed with respect to the K–L vectors is written in terms of modal coordinates as

η¨ l (t ) + 2ζl ωl η˙ l (t ) + ωl2 ηl (t ) = {φl }T {b}a(t )0 + {φl }T [Q ]{z (t )}0 η¨ lj (t ) + 2ζl ωl η˙ lj (t ) + ωl2 ηlj (t ) = {φl }T {b}a(t )j + {φl }T [Q ]{z (t )}j j = 1, . . . nkl , l = 1, . . . , nsm

(22)

where

{u(t )}0 =

nsm X

{φl }ηl (t )

{u(t )}j =

l =1

nsm X {φl }ηlj (t )

(23)

l=1 j

and ηl (t ) and ηl (t ) are modal coordinates, {φl } is the eigenvector associated with the eigenvalue ωl of the corresponding eigenvalue–eigenvector problem, ζl is a suitable damping ratio, and nsm is the number of structural modes considered in the dynamic analysis. The contribution of the remaining modes to the structural response can be approximated by using the static solution. This approach is reasonable due to the fact that the high frequency modes react essentially in a static manner when excited by low frequencies [27]. The contribution of the static solution can be incorporated directly in the present formulation. For simplicity in notation, however, such contribution is not considered in the presentation of the methodology. The equations for the temporal evolution of the mean and K–L vectors of the response process are solved by the second-order step by step integration scheme presented in Appendix A.

nkl X {u(t )}j ξj ,

8.4. Stochastic response sensitivity

j =1

It is seen from Eq. (12) that the approximations of the first two statistical moments of the stochastic performance functions involved in the compromise optimization problem are given in terms of their sensitivities. At the same time it is noted from Eqs. (6), (7) and (8) that the tradeoff analysis requires the same type of sensitivity information (gradients of the active constraints). In order to obtain such information a sensitivity analysis of the first two statistical moments of the stochastic performance functions is carried out. To introduce such analysis it is observed that the sensitivity of a stochastic response can be expressed with respect to the sensitivity of its K–L representation. For example, if σus (t , {θ }) represents the standard deviation of the displacement response process us (t ) (s component of the displacement vector {u(t )}), then

nkl X {˙u(t )} ≈ {˙u(t )}0 + {˙u(t )}j ξj j =1 nkl

{z (t )} ≈ {z (t )}0 +

8.3. Calculation of mean and K–L vectors

X {z (t )}j ξj

(20)

j =1

where {u(t )}0 , {˙u(t )}0 and {z (t )}0 denote the mean vectors, and {u(t )}j , {˙u(t )}j and {z (t )}j the jth K–L vectors of the response process, respectively. Once the Karhunen–Loève representation of the response process has been obtained, the response statistics can be calculated directly. In fact, the vectors {u(t )}0 , {˙u(t )}0 , {z (t )}0 , {u(t )}j , {˙u(t )}j and {z (t )}j uniquely specify the first two statistical moments of the linearized stochastic responses. For example, the second moment representation of the displacement vector {u(t )} can be calculated according to

{µ{u(t )} } = {u(t )}0 ,

nkl X

[ujs (t )]2 ,

(24)

j =1

and therefore nkl

[Γ{u(t )}{u(t )} (tn , tm )] =

σu2s (t , {θ }) =

X

{u(tn )}j {u(tm )}j

T

(21)

j =1

where {µ{u(t )} } is the mean vector, and [Γ{u(t )}{u(t )} ] is the covariance matrix of the displacement response process.

nkl P

nkl P

j

∂ u (t )

j

us (t ) ∂θs l

∂σus (t , {θ }) j =1 = s nkl ∂θl P

= j us

[ (t )]

j =1

2

j =1

j

j

∂ u (t )

us (t ) ∂θs l

σus (t , {θ})

,

(25)

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

591

Standard deviation

where

Peak time

design: {θ} + {Δθ}

{∂lu (t )}0 =

∂{u(t )}0 , ∂θl

{∂lu (t )}j =

∂{u(t )}j , ∂θl

{∂lz (t )}0 = {∂lz (t )}j =

∂{z (t )}0 ∂θl

∂{z (t )}j , ∂θl

(31)

On the other hand, the sensitivity of the linearized equations are given by design: {θ}

i ∂ h  E {h({u(t )}, {˙u(t )}, {z (t )})} (32) ∂θl {∂˙lz (t )}j = [Υ1 (t )]{∂lu (t )}j + [Υ2 (t )]{∂˙lu (t )}j + [Υ3 (t )]{∂lz (t )}j {∂˙lz (t )}0 =

Time

t*

nkl P

j

∂ u (t ∗ )

j

s us (t ∗ ) ∂θ l

∂ ∂σus (t , {θ}) j =1 max σus (t , {θ }) ≈ = s nkl ∂θl t ∂θl P ∗

[ujs (t ∗ )]2

nkl

P j =1

j

j

∂ u (t ∗ )

s us (t ∗ ) ∂θ l

(26)

(27)

To derive the equations for the sensitivity of the mean and K–L vectors, the K–L representation of the response process is substituted in Eqs. (16) and (19). Taking the derivative with respect to the design variable θl (l component of {θ}) yields

[M ]{∂¨ (t )} + [C ]{∂˙ (t )} + [K ]{∂ ( )} = {∂l (t )} [M ]{∂¨ (t )}j + [C ]{∂˙ (t )}j + [K ]{∂ ( )} = {∂l (t )}j , 0

u 0 l t u j l t

j = 1, . . . , nkl

(28)

∂[M ] ∂[C ] {¨u(t )}0 − {˙u(t )}0 ∂θl ∂θl

∂[K ] {u(t )}0 ∂θl ∂[M ] ∂[C ] {∂l (t )}j = [Q ]{∂lz (t )}j − {¨u(t )}j − {˙u(t )}j ∂θl ∂θl ∂[K ] − {u(t )}j , ∂θl −

η¨ ∂ l (t ) + 2ζl ωl η˙ ∂ l (t ) + ωl2 η∂ l (t ) = {φl }T {∂l (t )}0 j = 1, . . . nkl , l = 1, . . . , nsm

(35)

where η∂ l (t ) and η∂ l (t ) are modal coordinates, and all other terms have been previously defined. It is observed that the structure of the differential equations (32), (33) and (35) is similar to the equations that characterize the K–L representation of the response process (Eqs. (19) and (22)). Therefore, the iterative procedure proposed in Appendix A can also be used to integrate these equations. Finally, it is noted that the multi-objective optimization methodology presented in previous sections can also be applied to linear systems as a special case. The same procedure applies except that statistical equivalent linearization is not required in those cases. 9. Numerical examples 9.1. Example no. 1

0

with

{∂l (t )}0 = [Q ]{∂lz (t )}0 −

(34)

l =1

j

8.5. Sensitivity of mean and K–L vectors

u l u l

nsm X {φl }η∂j l (t )

η¨ ∂j l (t ) + 2ζl ωl η˙ ∂j l (t ) + ωl2 η∂j l (t ) = {φl }T {∂l (t )}j

In the above equation it has been assumed that the peak time is invariant in the neighborhood of a given design {θ}. Validation calculations have shown that such assumption is satisfied for the transient functions considered here, i.e. first and second-order statistics of the performance functions [25,28]. Fig. 3 shows a schematic illustration of this property. It is seen from the figure that even though the maximum values of the standard deviation for the designs {θ} and {θ} + {1θ }, where {1θ } represents a small perturbation, are different they are reached at the same time.

0

{∂lu (t )}j =

as

t

u l u l

nsm X {φl }η∂ l (t ), l =1

, σus (t ∗ , {θ}) where t ∗ is the peak time for σus (t , {θ }), that is   σus (t ∗ , {θ}) = max σus (t , {θ}) . =

where [∂l 1 (t )], [∂l 2 (t )] and [∂l 3 (t )] are the matrices that contain the sensitivity of the linearization coefficients. The sensitivity of the mass, damping and stiffness matrices required in Eqs. (29) and (30) can be evaluated directly, since these matrices can be generated at the element level in the same manner that the mass, damping and stiffness matrices are generated. In the above equations it has been assumed that the vector {b} and the matrix [Q ] are independent of the design variables, as stated before. To solve the equations for the sensitivity of the mean and K–L vectors, Eq. (28) is first written in terms of modal coordinates

{∂lu (t )}0 =

j =1

(33)

Υ

Υ

Υ

j

∂ u (t )

where ∂θs is the sensitivity of the sth component of the jth K–L l vector of the displacement response process. Similar expressions can be derived for the sensitivity of other stochastic responses. With this information the sensitivity of the maximum value of the stochastic response σus (t , {θ }) (required in Eq. (12)) is approximated as

Υ

Υ

Υ

+ [∂l 1 (t )]{u(t )}j + [∂l 2 (t )]{˙u(t )}j + [∂l 3 (t )]{z (t )}j

Fig. 3. Peak time invariant property for the statistics of the response functions.

The two-story frame structure with friction devices under earthquake loading shown in Fig. 4 is considered in the first example problem. The floor masses are mi = 1.3 × 106 kg, i = 1, 2, and initial linear interstory stiffnesses ki = 1.8 × 109 N/m, i = 1, 2. A 3% of critical damping is assumed in the model. For an improved earthquake resistance, the structure is reinforced with a nonlinear hysteretic device at each floor. It follows the inter-story restoring force law given by

(29)





r (t ) = kd δ(t ) − z 1 (t ) + z 2 (t ) ,

(30)

(36)

where kd denotes the initial stiffness of the device, δ(t ) is the relative displacement between floors, and z 1 (t ) and z 2 (t ) denote the plastic elongations of the device. Using the supplementary

592

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

a

r

b m kd

k2 2

k2 2

1

kd δ

m kd

k1 2

k1 2

a(t) Fig. 4. (a) Schematic sketch of two-story structure with nonlinear devices under ground acceleration. (b) Hysteresis loop between floors.

variable q(t ) = δ(t ) − z 1 (t ) + z 2 (t ), the plastic elongations are specified by the nonlinear differential equations [29]

˙ t ), q(t )), z˙ 1 (t ) = h1 (δ(

(37)

˙ t ), q(t )), z˙ (t ) = h2 (δ(

(38)

2

where



˙ t ), q(t )) = δ( ˙ t )H (δ( ˙ t )) H (q(t ) − qy ) h1 (δ(

q(t ) − qy qp − qy

H (qp − q(t ))

 + H (q(t ) − qp ) ,

(39)

and

−q(t ) − qy qp − qy  × H (qp + q(t )) + H (−q(t ) − qp ) , (40) 

˙ t ), q(t )) = −δ( ˙ t )H (−δ( ˙ t )) H (−q(t ) − qy ) h2 (δ(

where H (·) denotes the Heaviside step function, qy is a parameter specifying the onset of yielding, and kd qp is the maximum restoring force of the device. The values qp = 10−3 m, qy = 7.0 × 10−4 m, and kd = 6.0 × 108 N/m are used for each nonlinear element. The definition and evaluation of the corresponding linearization coefficients for this non-linear restoring force law are presented in Appendix B. At the same time, the basic relations for the evaluation of the sensitivity of the linearization coefficients are summarized in Appendix C. The earthquake induced ground acceleration is modelled as a non-stationary filtered white noise process. The filter is defined by the first-order differential equation 0 y1 (t ) d y2 (t ) −Ω12  = 0 dt y3 (t ) y4 (t ) Ω12





1



−2ξ1 Ω1 0

0 0 0

0 0 1

  

2ξ1 Ω1 −Ω22 −2ξ2 Ω2    y1 (t ) 0 y (t ) w(t )e(t ) ×  2 + , y3 (t ) 0 y4 (t ) 0



(41)

(42)

where w(t ) denotes white noise, and e(t ) denotes the envelope function defined as e(t ) =

e−0.5t − e−t max(e−0.5t − e−t )

,

0 < t < 10 s.

{F ({θ })} = {σ ({θ }), C ({θ })}.

(44)

For computational purposes the objective functions are normalized by their values at the initial design. The constraints for a given design are given by

β × σδi ({θ}) ≤ 0.002 m i = 1, 2 (45) where σδi ({θ }) is the maximum value of the standard deviation of the relative displacement between the (i − 1, i)th floors. The factor β may be interpreted as the safety index of the constraint condition. A larger factor requires the constraint to be more robust to the scatter of the stochastic process. A factor β = 3.0 is considered in this case. The optimization variables are the linear interstory stiffnesses k1 and k2 , with initial design k0i = 1.8 × 109 N/m, i = 1, 2, and side constraints 0.2 ≤ ki /k0i ≤ 2.0. The purpose of this example problem is to evaluate the performance of compromise programming in obtaining Pareto designs that best represent given aspiration levels. The corresponding compromise programming problem to the multi-objective optimization problem (44) and (45) can be written as Min d

and the ground acceleration is defined as a(t ) = Ω12 y1 (t ) + 2ξ1 Ω1 y2 (t ) − Ω22 y3 (t ) − 2ξ2 Ω2 y4 (t ),

The values Ω1 = 15.6 rad/s, ξ1 = 0.6, Ω2 = 1.0 rad/s, and ξ2 = 0.9, and white noise intensity I = 0.022 m2 /s3 have been used in this example. The sampling interval and the duration of the excitation are taken as 1t = 0.01 s, and T = 10 s, respectively. Then, the discrete-time white noise sequence ω(tj ) = √ I /1t ξj , (where ξj , j = 1, . . . , 1001, are independent, identically distributed standard Gaussian random variables) is considered in this case. The K–L representation of the filtered white noise process a(t ) can be obtained directly from the characterization of the filter response (see Appendix D). Due to the simplicity of the structural model all 1000 K–L vectors of the filtered white noise representation are used in this example. The multi-objective optimization problem is defined in terms of a cost function C ({θ }) and the maximum value of the standard deviation of the total displacement at the second floor of the structure σ ({θ }). For illustration purposes, the cost function is assumed to be proportional to the interstory stiffness of the first and second floor. Therefore the vector of objective functions is defined as

(43)

subject to C (k1 , k2 ) ≤ d(C (as) − C (id) ) + C (id)

σ (k1 , k2 ) ≤ d(σ (as) − σ (id) ) + σ (id) β × σδi (k1 , k2 ) ≤ 0.002 m i = 1, 2 0.2 ≤ ki /k0i ≤ 2.0,

d≥0

(46)

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

Fig. 5. Pareto points in the objective function space obtained by compromise programming.

Fig. 6. Iteration history of the compromise optimization problem in terms of the auxiliary variable d (example 1).

where the small parameter α (see Eq. (2)) has been taken as zero in this case. In the above equations σ (as) and C (as) represent the aspiration level of the objective functions σ and C , respectively, while σ (id) and C (id) are the corresponding ideal values (individual optima for each of the objective functions). The ideal point in the objective function space (σ (id) , C (id) ) = (1.36 × 10−3 m, 1.76 × 109 N/m) as well as the points corresponding to single objectives, i.e, minimum response design problem (σ , C ) = (1.36 × 10−3 m, 7.20 × 109 N/m), and minimum cost design problem (σ , C ) = (3.96 × 10−3 m, 1.76 × 109 N/m) are shown in Fig. 5. For illustration purposes two aspiration levels are considered

(σ (as) , C (as) )1 = (2.20 × 10−3 m, 2.00 × 109 N/m) (σ

(as)

,C

(as)

)2 = (3.30 × 10

−3

m, 6.00 × 10 N/m). 9

(as)

Fig. 7. Iteration history of the compromise optimization problem in terms of the objective function σ (example 1).

Fig. 8. Iteration history of the compromise optimization problem in terms of the objective function C (example 1).

objective function σ , and the objective function C is shown in Figs. 6–8, respectively. The compromise programming problem associated with the aspiration level (σ (as) , C (as) )1 is considered in these figures. The initial design (initial point) for the compromise programming problem is taken equal to the aspiration level (σ (as) , C (as) )1 . It is seen that the process converges in less than four iterations. Recall that each iteration requires just one stochastic analysis (sensitivity analysis). Therefore, the Pareto point is obtained in a very efficient manner. The same efficiency is found for other aspiration levels. Thus, the method is robust in terms of the aspiration levels considered in the analysis. 9.2. Example no. 2

(47)

(as)

Note that the aspiration level (σ , C )1 is unfeasible while (σ (as) , C (as) )2 is feasible. As previously pointed out, these aspiration levels are defined by the designer prior to solving the compromise programming problem. By solving the compromise optimization problem the following designs are obtained

(σ , C )1 = (3.35 × 10−3 m, 2.33 × 109 N/m) (σ , C )2 = (2.23 × 10−3 m, 4.00 × 109 N/m).

593

(48)

For comparison, the Pareto front is also shown in the figure. It is seen that the points (σ , C )1 and (σ , C )2 are Pareto points in the sense that no other solutions are superior to them when the two objective functions are considered simultaneously. It is also observed that the Pareto points minimize the distance between the ideal point and the Pareto front in the direction defined by the ideal point and the aspiration levels. The iteration history of the optimization process corresponding to the compromise optimization problem in terms of the auxiliary variable d, the

9.2.1. Description A finite element building model is considered in the second numerical example. The elevation view of the model is shown in Fig. 9. Each floor is modeled by rectangular plate elements made of reinforced concrete with a thickness of t = 0.18 m. Material properties of the reinforced concrete have been assumed as follows: Young’s modulus E = 2.1 × 1010 N/m2 , Poisson ratio µ = 0.25, and weight density ρ = 2.5 ton/m3 . The total mass of each floor is given by 8.8 × 104 kg. All floors have a constant height equal to 3.2 m, leading to a total height of 9.6 m. Additionally, twenty four steel brace elements are used in the model, with material properties E = 2.0 × 1011 N/m2 , and weight density ρ = 7.42 ton/m3 . A 3% of critical damping for the first modes is considered in the model. For aseismic design purposes the structure is reinforced with non-linear hysteretic devices (shear panels) as shown in the figure. These elements provide additional resistance against relative displacements between floors. The devices follow the interstory restoring force law r (t ) = k(1(t ) −

594

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599 VARIABLE 6

VARIABLE 3

VARIABLE 6 VARIABLE 3

VARIABLE 2

VARIABLE 2

SHEAR PANEL

VARIABLE 5

VARIABLE 1

SHEAR PANEL

10.9 m SHEAR PANEL

VARIABLE 1

7.5 m

7.5 m VARIABLE 4 10.9 m

6.0 m

z y x

Fig. 9. Finite element building model.

z 1 (t ) + z 2 (t )), where k denotes the initial stiffness of the nonlinear element, 1(t ) is the relative displacement between the degrees of freedom where the element is located, and z 1 (t ) and z 2 (t ) denote the plastic deformation of the non-linear element. The restoring force r (t ) acts with the same orientation as the relative displacement 1(t ). The evolution of the plastic deformations z 1 (t ) and z 2 (t ) with respect to time is given by the non-linear differential Eqs. (37)–(40). Because of the yielding, energy dissipation due to hysteresis is introduced in the structural response. The devices have initial stiffness k = 2.4 × 108 N/m, and model parameters qp = 0.0028 m, and qy = 0.00194 m. The structure is excited horizontally by earthquake excitation a(t ) in the x direction. It is modeled as in the previous example with model parameters Ω1 = 15.6 rad/s, ξ1 = 0.6, Ω2 = 1.0 rad/s, ξ2 = 0.995, and white noise intensity I = 0.36 m2 /s3 , with time modulation function e(t ) defined as e(t ) = e−0.25t − e−0.5t /max(e−0.25t − e−0.5t ), 0 < t < 20 s. The K–L representation of the filtered white noise can be obtained as before, that is, according to the method presented in Appendix D. The process a(t ) is represented by the K–L expansion using a time span of [0, 20] s with a time step 1t = 0.01 s, yielding a series expansion with a total of 2001 terms. The actual number of K–L vectors used in this case is equal to nkl = 400. On the other hand, 12 out of the total number of classical modes (considering the static contribution of the higher modes) are retained in the analysis. These numbers proved out to be adequate for the example under consideration.

functions while satisfying design constraints given in terms of the interstory drift ratios at the center of the building. The multiobjective optimization problem is written as

9.2.2. Multi-objective optimization problem The multi-objective optimization problem includes the following objective functions: the structural weight represented by the weight of the brace elements (W ), and the maximum value of the standard deviation of the total displacement at the center of the third floor (σ ). The design variables of the optimization problem are the areas of the brace elements. Each floor is linked to two design variables which are related to resistant planes in the x and y direction, respectively (see Fig. 9). Thus, the vector of design variables {θ } has six components with initial design θ1 = 24.0 cm2 , θ2 = 20.0 cm2 , θ3 = 12.0 cm2 , θ4 = 24.0 cm2 , θ5 = 20.0 cm2 , and θ6 = 12.0 cm2 , and side constraints 15.0 cm2 ≤ θ1 , θ4 ≤ 45.0 cm2 , 10.0 cm2 ≤ θ2 , θ5 ≤ 35.0 cm2 , and 6.0 cm2 ≤ θ3 , θ6 ≤ 25.0 cm2 . The goal of the design problem is to find the size of the brace elements that achieves a compromise among the two objective

subject to

Minimize {F ({θ})} = {σ ({θ }), W ({θ})} subject to

β × σ1i ≤ σ ∗ ,

i = 1, 2, 3

(49)

where σ1i is the maximum value of the standard deviation of the relative displacement between the (i − 1, i)th floors at the center of the building, and σ ∗ = 0.016 m is the upper bound for the interstory drift ratio constraints (0.5% of the interstory height). A factor β = 2.5 (safety index) is considered in this case. The structural weight is given explicitly in terms of the design variables as W ({θ}) = 20.2θ1 + 20.2θ2 + 20.2θ3 + 33.7θ4 + 33.7θ5 + 33.7θ6 . 9.2.3. Pareto point The points corresponding to single objective problems, i.e., minimum response design problem Minimize σ ({θ}) subject to

β × σ1i ≤ σ ∗ ,

i = 1, 2, 3

(50)

and minimum structural weight design problem Minimize W ({θ})

β × σ1i ≤ σ ∗ ,

i = 1, 2, 3

(51)

as well as the ideal point, that is, the vector that contains the individual optima for each of the objective functions (σ (id) = 0.00541 m, W (id) = 1822 kg) are shown in Fig. 10. In this figure the values of the standard deviation of the total displacement at the third floor σ are normalized by the height of the floors multiplied by 10−3 , i.e., σˆ = σ /(3.2 × 10−3 ). On the other hand, the initial aspiration levels are assumed to be W (as) = 1.03 × W (id) , and σ (as) = 1.03 × σ (id) . Note that these aspiration levels are very close to the ideal solution and therefore it is expected that this design is unfeasible. In fact, it can be shown that the design corresponding to these aspiration levels does not satisfy the constraints of the multi-objective optimization problem (49). By

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

Fig. 10. Pareto point in the objective function space obtained by compromise programming.

595

Fig. 12. Trajectory of the optimizer for the compromise optimization problem in the objective function space.

Table 1 Compromise design. Design variables

Pareto design

θ1 cm2 θ2 cm2 θ3 cm2 θ4 cm2 θ5 cm2 θ6 cm2

19.3 14.2 10.7 15.0 10.0 6.0 3

No. of stochastic analyses

Table 2 Design constraints. Fig. 11. Iteration history of the compromise optimization problem in terms of the auxiliary variable d (example 2).

solving the corresponding compromise programming problem a Pareto point that best represents the designer initial preferences is obtained. The optimum objective function values are σ ∗ = 0.00573 m, and W ∗ = 1938 kg. As previously pointed out this point belong to the set of solutions which are optimal in the sense that no other solutions are superior to them when the two objective functions are considered simultaneously. The corresponding design for this Pareto point is given in Table 1. The values of the constraints at the final design are shown in Table 2. It can be seen that the constraints related to the interstory drift ratios of the first and second floor are active at the optimum. At the Pareto point the design variables θ4 , θ5 , and θ6 , are equal to their lower bounds. This is reasonable since the excitation is acting in the x direction and the model is symmetric about that axis. Therefore, the resistant planes in the y direction do not work under the stochastic load. From Figs. 11 and 12 it is seen that the optimization process of the compromise programming problem converges in less than three iterations. These figures show the iteration history of the optimization process in terms of the auxiliary variable d (objective function of the compromise optimization problem), and in terms of the objective functions in the objective function space, respectively. The initial design (initial point) for the compromise optimization problem is taken equal to (W = 2177 kg, σˆ = 1.72). The time needed for finding the Pareto solution of this highly demanding problem is less than one hour. Therefore the compromise design that best represents the designer preferences is obtained in a reasonable time. Additional validation calculations have shown that the efficiency of the method remains invariant regardless of the selected aspiration levels.

Constraints (m)

Pareto design

σ11 σ12 σ13

0.016a 0.016a 0.011

a

Active constraints.

is normalized as before, that is, σˆ = σ /(3.2 × 10−3 ). First, the feasible directions with the greatest improvement of the objective functions σˆ and W are given by (see Eq. (6))

{dσˆ } = −[P ]∇ σˆ {θ ∗ } {dW } = −[P ]∇ W ({θ ∗ }) where

[P ] = [I ] − [B]T ([B][B]T )−1 [B]

(53)

is the projection matrix. The rows of the matrix [B] are defined by the gradients of the active constraints including the active bounds, in this case

 ∇ σˆ 11 ({θ ∗ })T ∇ σˆ 1 ({θ ∗ })T  2   ∗ T  [B] =   ∇ g4 ({θ ∗ })T   ∇ g5 ({θ })  ∇ g6 ({θ ∗ })T 

(54)

where σˆ 11 and σˆ 12 are the normalized standard deviation of the relative displacements between floors, and the constraints g4 ≤ 0, g5 ≤ 0 and g6 ≤ 0 represent the active bound constraints, with g4 ({θ ∗ }) = −θ4 + 15,

9.2.4. Sensitivity analysis The Pareto sensitivity analysis is conducted at the current design {θ ∗ } to study the effect of changes in one objective function on the other objective. In what follows, the objective function σ

(52)

g6 ({θ }) = −θ6 + 6. ∗

g5 ({θ ∗ }) = −θ5 + 10,

(55)

Recall that the design variables θ4 , θ5 , and θ6 , are equal to their lower bounds at the Pareto point (see Table 1). Numerically the

596

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

matrix [B] is equal to

−0.101  −0.02 [B] =  0.00 0.00 0.00

−0.027 −0.179 0.00 0.00 0.00

−0.037 −0.036 0.00 0.00 0.00

0.00 0.00 −1.00 0.00 0.00

0.00 0.00 0.00 −1.00 0.00

0.00  0.00  0.00  . 0.00 −1.00

(56) From the coefficients of this matrix it is clear that the area of the brace elements of the first floor controls the response of the interstory drift ratio of the first floor, as expected. Similarly, the area of the brace elements of the second floor dominates the response of the interstory drift ratio of the second floor. The negative sign of the coefficients of the first two rows indicates that an increase in the area of the brace elements is related to a decrease in the interstory drift ratios. The feasible directions {dσˆ } and {dW } are then obtained as

−0.023

 2.96   1.51  −9.28  {dW } =   0 .  

−0.012  0.072   {dσˆ } =   0 ,   0 0

(57)

0 0

These directions correspond to the projection of the steepest descent direction of the objective functions onto the tangent plane defined by the active constraints at the Pareto point (the constraints related to the interstory drift ratios of the first and second floor, and the side constraints for the design variables θ4 , θ5 and θ6 ). Note that the last three components of these vectors are zero since the corresponding design variables are fixed at their lower bound values and therefore they can not change. Accordingly, the values of the objective functions can change only with respect to θ1 , θ2 and θ3 . The corresponding sensitivity information along the feasible descent direction of each of the objective functions is given by the tradeoff matrix

  [T ] =

1 δ W /δ σˆ

 =

1 −128.2

δ σˆ /δ W 1



1

 = 

 −0.0078 1

.

{dW }T {dσˆ } {dσˆ }T {dσˆ }

 {dσˆ }T {dW } {dW }T {dW }    1

(58)

It can be seen that an increase in the structural weight of the system is associated with a decrease in the maximum value of the standard deviation of the total displacement at the top of the building. Similarly, an increase in the displacement response is achieved by a decrease in the structural weight of the model. The fact that the off-diagonal terms of the tradeoff matrix are negatives indicates that the current design is Pareto optimum. 9.2.5. Tradeoff analysis Fig. 13 shows the first-order approximation to the Pareto surface around the present design {θ ∗ }. The approximate surface allows to explore in an effective manner the Pareto surface locally. For example, new aspiration levels can be generated directly from the approximation, which can be used as approximate new Pareto points. The new aspiration levels can also be used to obtain non-inferior solutions by solving the corresponding compromise optimization problems. Fig. 13 also shows the accuracy of the local Pareto surface approximation by comparing the predicted points with the exact ones in the vicinity of the Pareto point. The exact solution is obtained by solving the compromise optimization

Fig. 13. Comparison between Pareto surface approximation and exact Pareto points in the vicinity of the computed Pareto point.

problem associated with the predicted point. It is seen that the values of the objective functions are very similar. Therefore, the local Pareto surface approximation is sufficiently accurate in this case to be used as an efficient exploration tool around the Pareto point. Then, the aspirations predicted by the approximation are attainable and approximately on the Pareto front. This means that the new approximate Pareto points are closed to the designer preferences. This is important for the rapid convergence of the decision-making process and for the confidence of the analyst. If the designer is not satisfied with the current set of compromise designs, new aspiration levels and the corresponding approximate Pareto surface can be obtained. In this manner the designer can conduct a decision-making analysis by means of an efficient interactive procedure. 10. Conclusions The goal of the present work was to introduce an efficient computational procedure for multi-objective optimization of non-linear structural systems under stochastic excitation. For that purpose several techniques such as equivalent statistical linearization, orthogonal expansion of stochastic processes, modal analysis, sensitivity analysis, and mixed linearization have been integrated. Numerical results show that the total number of stochastic analyses required during the entire design process is in general very small. Hence, different compromise solutions including the design that best represents the outcome that the designer considers potentially satisfactory are obtained in an efficient manner. At the same time an effective sensitivity analysis of Pareto solutions can be carried out by the proposed formulation. Such sensitivity information provides the designer a practical tool for efficient exploration around Pareto solutions and for decision making and tradeoff analysis. The results indicate that the Pareto designs obtained by the Pareto surface approximation approach are closed to the exact ones for the cases considered here. Using this approach the aspiration levels that reflect the preferences of the designer can be easily generated, and the actual Pareto point obtained from these aspirations is close to the aspirations themselves. This is important from a practical point of view since the designer can obtain a satisfactory final Pareto design in a relatively small number of design iterations. It is concluded that the proposed multi-objective optimization procedure can be very effective in realistic engineering problems such as non-linear structural systems under stochastic excitation. Acknowledgment The research reported here was supported in part by CONICYT under grant number 1070903 which is gratefully acknowledged by the author.

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

597

Appendix A. Numerical integration for the mean and K–L vectors

where E (·) denotes the expectation operator, and all other terms have been previously defined. The linearization coefficients can be determined by the integrals

Provided that the response process at time t is known, the following iterative procedure is implemented to compute the response at time t + 1t. At the beginning of the iteration (s = 0) it is assumed that

Z qp Z ∞   q − qy ˙ q)dδ˙ dq ˙ q) = f (δ, λz1 ,0 = E h1 (δ, δ˙ ˙ 0 qp − qy q=qy δ= Z ∞ Z ∞ ˙ q)dδ˙ dq δ˙ f (δ, +

{z (t + 1t )}0(s=0) = {z (t )}0 ,

{˙z (t + 1t )}0(s=0) = {˙z (t )}0

{z (t + 1t )}j(s=0) = {z (t )}j ,

{˙z (t + 1t )}j(s=0) = {˙z (t )}j . (A.1)

Then, the decoupled equations (22) are integrated by the (s) Newmark constant acceleration method, to obtain ηl (t + 1t ), j(s) l

and η (t + 1t ), l = 1, . . . , nsm . With this information, Eq. (19) is solved by a second-order step by step integration scheme

{z (t + 1t )}0(s) = {z (t )}0 + {z (t + 1t )}j(s) = {z (t )}j +

1t

2 1t 2

{˙z (t )}0 + {˙z (t + 1t )}0(s−1) {˙z (t )}j + {˙z (t + 1t )}


j(s−1)

j = 1, . . . , nkl ,

,

λz1 ,2

 Z zp Z ∞ ˙ q) q − qy ∂ h1 (δ, = f (δ˙ , q)dδ˙ dq =E q ˙ ∂ δ˙ q=qy δ=0 p − qy Z ∞ Z ∞ ˙ q)dδ˙ dq f (δ, +

˙ q) = f (δ,

(A.3)

(A.4)

where the matrices [Υ1 (t + 1t )]j(s) , [Υ2 (t + 1t )]j(s) and [Υ3 (t + 1t )]j(s) contain the linearization coefficients computed from the current solution at the sth iteration within the time interval [t , t + 1t ], that is,

{u(t + 1t )}0(s) , {˙u(t + 1t )}0(s) , {z (t + 1t )}0(s−1) j(s)

j(s)

{u(t + 1t )} , {˙u(t + 1t )} , {z (t + 1t )}

j(s−1)

.

(A.5) (A.6)

Similarly, the mean value of the vector function h, i.e., E ({h({u(t + 1t )}, {˙u(t + 1t )}, {z (t + 1t )})})

k{z (t + 1t )} − {z (t + 1t )} k{z (t + 1t )}j(s−1) k

j(s−1)

k

1 2π σq σδ˙

q

− 2ρδ, ˙q

2 1 − ρδ, ˙q



q − µq

exp −



σq

"

1

q − µq

σq

2 2(1 − ρδ, ˙ q)

δ˙ − µδ˙ σδ˙



 +

2

δ˙ − µδ˙ σδ˙

2 #

, (B.5)

where µq and µδ˙ denote the mean and σq and σδ˙ denote the standard deviation of the random variables q and δ˙ , respectively, and ρδ, ˙ q denotes the correlation coefficient. These statistics can be computed in terms of the mean and K–L vectors of the responses {u(t )}, {˙u(t )} and {z (t )}, since the relative displacement δ(t ) can be expressed with respect to the components of the displacement vector {u(t )}. Some of the terms involved in the above integrals can be integrated analytically. The remaining terms can be expressed with respect to the error function, and therefore they can be computed efficiently by using a suitable numerical integration scheme. The statistical linearization of the equation z˙2 (t ) = h2 (δ˙ (t ), q(t )) can be carried out accordingly.

(A.7)

is computed from the current solution of the response. The procedure is repeated until j(s)

(B.4)

˙ 0 δ=

˙ q) has been assumed to be jointly where the response process (δ, Gaussian distributed with probability density function

0(s−1)

j = 1, . . . , nkl ,

(B.3)



q=qp

where the quantities {z (t )}0 , {˙z (t )}0 , {z (t )}j and {˙z (t )}j are known from the previous time step (solution at time t), and

+ [Υ3 (t + 1t )]j(s) {z (t + 1t )}j(s−1) ,

 Z qp Z ∞ ˙ q) ∂ h1 (δ, δ˙ ˙ q)dδ˙ dq f (δ, = ∂q ˙ 0 qp − qy q=qy δ=

and




(A.2)

{˙z (t + 1t )} = E ({h({u(t + 1t )}, {˙u(t + 1t )}, {z (t + 1t )})}) {˙z (t + 1t )}j(s−1) = [Υ1 (t + 1t )]j(s) {u(t + 1t )}j(s) + [Υ2 (t + 1t )]j(s) {˙u(t + 1t )}j(s)



λz1 ,1 = E

(B.2)

˙ 0 δ=

q=qp

≤  j = 0, 1, . . . , nkl (A.8)

where  is a user-defined small scalar. Numerical experience has shown that in general only few iterations (less than five) are required within each time interval [t , t + 1t ] if small time steps, compared with the periods of the modal equations, are used in the integration [18]. The relationships given by Eqs. (A.3) and (A.4) are evaluated in local component specific coordinates increasing in this manner the efficiency of the integration scheme for calculating the K–L representation of the response process.

Appendix C. Sensitivity of linearization coefficients Using the definition of the linearization coefficients given in Appendix B, the sensitivity of the coefficients can be determined as follow

∂ λz1 ,0 =

∂λz1 ,0 = ∂θl Z ∞ Z +

∂

∂λz1 ,1 = = ∂θl

Z

δ˙

qp

∞

Z

q=qy

δ˙

˙ q) q − qy ∂ f (δ, qp − qy

∂θl

dδ˙ dq

˙ q) ∂ f (δ, dδ˙ dq ∂θl

∞

Z

∞ ˙ 0 δ=

q =q y

˙ 0 δ=

q =q p

λz1 ,1

qp

Z

˙ 0 δ=

(C.1)

˙ q) δ˙ ∂ f (δ, dδ˙ dq qp − qy ∂θl

(C.2)

and Appendix B. Linearization coefficients Since the non-linear differential equations (39) and (40) are of the same structures it suffices to consider, for example, the ˙ t ), q(t )). The linearization of the non-linear equation z˙1 (t ) = h1 (δ( equation is replaced by the linear relation z˙1 (t ) = λz1 ,0 + λz1 ,1 [q(t ) − E (q(t ))] + λz1 ,2 [δ˙ (t ) − E (δ˙ (t ))], (B.1)

∂ λz1 ,2 =

∂λz1 ,2 = ∂θl Z ∞ Z + q =q p

Z

qp q =q y

∞ ˙ 0 δ=

Z

∞ ˙ 0 δ=

˙ q) q − qy ∂ f (δ, qp − qy

˙ q) ∂ f (δ, dδ˙ dq ∂θl

∂θl

dδ˙ dq

(C.3)

˙ q) has been assumed to where as before the response process (δ, be jointly Gaussian distributed. Given the algebraic structure of the

598

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599

˙ q) (see (B.5)) it can be shown that probability density function f (δ,

with initial condition

˙ q) ∂ f (δ, ˙ q) f (δ, ˙ q) = P2 (δ, ∂θl

[Γ{y}{y} (tn = tm , tm )] = [Γ{y}{y} (tm )].

(C.4)

˙ q) is a multivariable polynomial function of secondwhere P2 (δ, order defined as

˙ q) = τ0 + τ1 (q − µq ) + τ2 (δ˙ − µδ˙ ) + τ3 (q − µq )(δ˙ − µδ˙ ) P2 (δ, + τ4 (q − µq )2 + τ5 (δ˙ − µδ˙ )2

(C.5)

and where the polynomial coefficients are given by

τ0 = τ1 = τ2 = τ3 = τ4 =

−1 q

2π σq σδ˙

2 1 − ρδ, ˙q

 q ∂  2 2π σq σδ˙ 1 − ρδ, ˙q ∂θl

ρδ, ∂   ∂   ˙q − µ µδ˙ (C.7) q 2 2 σq2 (1 − ρδ, σq σδ˙ (1 − ρδ, ˙ q ) ∂θl ˙ q ) ∂θl ρδ, 1 ∂   ∂   ˙q (C.8) µδ˙ − µq 2 2 2 σδ˙ (1 − ρδ, σq σδ˙ (1 − ρδ, ˙ q ) ∂θl ˙ q ) ∂θl ! ρδ, ∂ ˙q , 2 ∂θl σq σδ˙ (1 − ρδ, ˙ q) ! ! (C.9) ∂ −1 −1 ∂ , τ5 = . 2 2 ∂θl 2σq2 (1 − ρδ, ∂θl 2σδ˙2 (1 − ρδ, ˙ q) ˙ q) 1

The partial derivatives involved in the definition of the polynomial coefficients can be expressed directly in terms of the sensitivity of the mean and K–L vectors of the responses {u(t )}, {˙u(t )} and {z (t )}, since as previously pointed out, the relative displacement δ(t ) can be expressed with respect to the components of the displacement vector {u(t )}. As in the case of the linearization coefficients (Appendix B) some of the terms involved in the above integrals can be integrated analytically. The remaining terms can be expressed in terms of the error function, and therefore they can be evaluated directly by using an appropriate numerical integration scheme. Appendix D. K–L representation of filtered white noise The covariance matrix [Γ{y}{y} (t )] describing the correlation of the filter parameters for each point in time t, is described by the Lyapunov matrix differential equation [30]

[Γ˙ {y}{y} (t )] = [A][Γ{y}{y} (t )] + [Γ{y}{y} (t )][A]T + {b(t )}I {b(t )}T ,

(D.1)

with initial condition

[Γ{y}{y} (t = 0)] = [0]

(D.2)

where



0

− [A] = 

Ω12

1

−2ξ1 Ω1

0

0

Ω12

2ξ1 Ω1

0 0 0

−

Ω22

0 0 1

  ,

(D.3)

−2ξ2 Ω2

{b(t )}T

= {0, e(t ), 0, 0}, and I is the white noise intensity. In addition, the covariance matrix [Γ{y}{y} (tn , tm )] describing the correlation of the filter parameters between the two time instances tn and tm can be calculated according to the matrix differential equation [30] ∂ [Γ{y}{y} (tn , tm )] = [A][Γ{y}{y} (tn , tm )] ∂ tn

From the definition of the stochastic process a(t ) in terms of the state-space variables of the filter {y(t )} (see Eq. (42)), the components of the covariance matrix of the discrete process [Γ{a}{a} ] can be calculated directly. For example E (a(tn )a(tm )) = {γ }T [Γ{y}{y} (tn , tm )]{γ }

(D.4)

(D.6)

where

{γ }T = {Ω12 , 2ξ1 Ω1 , −Ω22 , −2ξ2 Ω2 }. (C.6)

(D.5)

(D.7)

By solving the algebraic eigenproblem

[Γ{a}{a} ]{ψj } = λj {ψj }

(D.8)

of the covariance matrix [Γ{a}{a} ], the K–L representation of the discrete filtered white noise process {a} is determined as

{a} ≈ {a}0 +

nkl X {a}j ξj

(D.9)

j =1

where

{a}T = {a(t0 ), a(t1 ), . . . ., a(tN )},

(D.10)

N is the number of time points, nkl is the order of truncation of the series expansion, {a}0 is the mean vector of the process, and {a}j , j = 1, . . . ., nkl are the K–L vectors defined as

{ a} j =

p

λj {ψj },

j = 1, . . . , nkl

(D.11)

References [1] Steuer RE. Multiple criteria optimization: Theory, computation and applications. New York: Wiley; 1986. [2] Pareto V. Cour d’Economie Politique. Geneva: Librarie Droz; 1964. [3] Koski J. Multicriteria truss optimization. In: Stadler W, editor. Multicriteria optimization in engineering and sciences. New York: Plenum; 1988. p. 263–306. [4] Stadler W, Daure J. Multicriteria optimzation in engineering: A tutorial and survey. In: Kamat MP, editor. Structural optimization: Status and promise. Progress in astronautics and aeronautics, vol. 150. Washington (DC): AIAA; 1992. p. 209–44. [5] Miettinen KM. Nonlinear multiobjective optimization. In: International series in operations research and management science. Norwell (MA): Kluwer Academic; 1999. p. 29–43. [6] Rao SS, Venkayya VB, Khot NS. Game theory approach for the integrated design of structures and controls. AIAA 1988;26(4):654–67. [7] Ulrich KT, Eppinger SD. Product design and development. 2nd ed. New York: McGraw-Hill; 2000. p. 138–59. [8] Papadrakakis M, Lagaros ND, Plevris V. Multiobjective optimization of space structures under static and seismic loadings conditions. Engineeering Optimization 2002;34:645–69. [9] Messac A, Mattson CA. Generating well-distributed sets of pareto points for engineering design using physical programming. Optimization and Engineering 2002;3:431–50. [10] Tappeta RV, Renaud JE. Interactive multiobjective optimization procedure. AIAA Journal 1999;37(7):881–9. [11] Goldberg DE. Genetic algorithms in search, optimization and machine learning. Reading (MA): Addison-Wesley; 1989. [12] Zitzler E, Thiele L. Multiobjective evolutionary algorithms: A comparative case study and the strength pareto approach. IEEE Transaction on Evolutionary Computation 1999;3(4):257–71. [13] Papadimitriou C. Pareto optimal sensor locations for structural identification. Computer Methods in Applied Mechanics and Engineering 2005;194: 1655–73. [14] Lagaros ND, Charmpis DC, Papadrakakis M. An adaptive neural network strategy for improving the computational performance of evolutionary structural optimization. Computer Methods in Applied Mechanics and Engineering 2005;194:3374–93. [15] Sawaragi Y, Nakayama H, Tanino T. Theory of multiobjective optimization. Orlando (FL): Academic; 1985. [16] Belegundu AD, Chandrupatla TR. Optimization concepts and applications in engineering. New Jersey: Prentice Hall; 1999. [17] Luenberger DG. Introduction to linear and nonlinear programming. AddisonWesley; 1973.

H.A. Jensen / Probabilistic Engineering Mechanics 24 (2009) 585–599 [18] Jensen HA. Structural optimization of non-linear systems under stochastic excitation. Journal of Probabilistic Engineering Mechanics 2006;21:397–409. [19] Fleury C, Braibant V. Structural optimization: A new dual method using mixed variables. International Journal for Numerical Methods in Engineering 1986; 23(3):409–28. [20] Haftka RT, Gürdal Z. Elements of structural optimization. 3rd ed. The Netherland: Kluwer Academic Publishers; 1992. [21] Schmit LA. Structural synthesis—Its genesis and development. AIAA Journal 1981;19(8):1249–63. [22] Loève M. Probability theory. 3rd ed. Princeton (New Jersey): D. Van Nostrand Company, Inc.; 1963. [23] Roberts JB, Spanos PD. Random vibration and statistical linearization. New York: Wiley; 1990. [24] Lutes LD, Sarkani S. Random vibrations. Oxford (UK): Elsevier; 2004.

599

[25] Jensen HA. Tradeoff analysis of dynamical systems under stochastic loading. In: Proceedings 7th european conference on structural dynamics EURODYN 2008 paper E33. 2008. [26] Emam HH, Pradlwater HJ, Schuëller GI. A computational procedure for the implementation of equivalent linearization in finite element analysis. Earthquake Engineering and Structural Dynamic 2000;29:1–17. [27] Bathe KJ. Finite element procedures. New Jersey: Prentice Hall; 1996. [28] Jensen HA, Marillanca A, Peñaloza O. A computational procedure for response statistics-based optimization of stochastic non-linear FE models. Computer Methods in Applied Mechanics and Engineering 2008;198:125–37. [29] Pradlwarter HJ. Deterministic integration algorithms for stochastic response computations of FE-systems. Computers & Structures 2002;80:1489–502. [30] Soong TT, Grigoriu M. Random vibration of mechanical and structural systems. Englewood Cliffs (NJ): Prentice-Hall, Inc; 1993.