Static and dynamic topology optimization: an innovative unifying approach

This paper presents a topology optimization approach that is innovative with respect to two distinct matters. First of all the proposed formulation is capable to handle static and dynamic topology optimization with virtually no modifications. Secondly, the approach is inherently a multi-input multi-output one, i.e., multiple objectives can be pursued in the presence of multiple loads. The input-to-output transfer matrix, say G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\varvec{G}}}}$$\end{document}, is the key ingredient that governs the algebraic mapping between applied loads and structural response. In statics G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\varvec{G}}}}$$\end{document} depends on the design variables only, whereas it depends on the frequency variable as well in the dynamic case. The Singular Value Decomposition (SVD) of G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\varvec{G}}}}$$\end{document} represents then the core of the proposed approach. Singular values are shown to be the gains of the input/output mapping and are used to compute proper norms of G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\varvec{G}}}$$\end{document} that represent the goal functions to be minimized. Singular vectors provide at no extra cost the plant directions, i.e., the load combination factors that stress the structure the most. Numerical examples are discussed in much detail and open issues object of ongoing investigations are highlighted. A full Matlab code handling the static topology optimization problem is provided as an online Appendix to the manuscript. Its extension to the dynamic case may be gathered following the formulation proposed in Sect. 5.


Introduction and motivation
Topology optimization is nowadays a mature discipline from both engineering and computational viewpoints. A systematic view about principles, computational methods, and applications was given in the classic book (Bendsøe and Sigmund 2003) that dates back some twenty years ago. The vast majority of contributions in both static and dynamic regimes aim at the minimization of (possibly slightly different definitions of) the structural compliance in the presence of further goals and/or constraints.
With no aim of completeness, as far as static response is concerned stress constraints were investigated in Bruggi and Duysinx (2012), Bruggi and Venini (2008), Giraldo-Londono et al. (2022), and Collet et al. (2017) focused on fatigue response, design of compliant mechanisms was tackled in Emmendoerfter et al. (2022), whereas Ferrari and Sigmund (2020b) and Gersborg-Hansen et al. (2006) addressed heat conduction and buckling, respectively. The development of fast computational methods was pursued in Andreassen et al. (2011) and Ferrari and Sigmund (2020a), whereas the advantages of mesh adaptivity were exploited in Salazar de Troya and Tortorelli (2018), just to mention a few. The above list is by no means complete and reference is made to the references therein and to the review paper (Sigmund and Maute 2013) for a comprehensive list of contributions and applications.
When it comes to dynamic response, a first main subdivision should be made between time-domain and frequencydomain approaches. As for time-domain approaches, Min et al. (1999) should be cited as a pioneering contribution, Giraldo-Londono and Paulino (2021) represents a milestone in that a general approach is presented along with a Matlab implementation that fits a more general modular environment, whereas (Jensen et al. 2014) faces and solves the fundamental issue of sensitivity computation. Switching to frequency-domain methods, Olhoff and Du (2016) presents an incremental approach for the minimization of the dynamic compliance in the presence of prescribed low or high value of the excitation frequency, Yang and Li (2014) minimizes the dynamic compliance of structures at resonance frequencies in thermal environments, whereas the frequent case of harmonic excitation is investigated in Liu et al. (2015).
Given the scenario above, this paper introduces an innovative formulation for the topology optimization of Multi-Input (i.e., multi-loaded) Multi-Output (i.e., multi-objective) structures. The proposed framework includes classic compliance optimization as a representative case but is general enough to allow extensions to alternative definitions of compliance as well as to the adoption of other merit functions that may be more appropriate for specific applications. The main idea behind the proposed approach is to write explicitly the input/output transfer matrix, say G , that maps the acting loads to the structural response. Using the adopted format, the structural response vector (that more often than not goes under the name of output vector) can be any linear combination of the state vector, i.e., it may include any nodal displacement of the discretized structure. Interestingly enough, the structural compliance (that is the goal function used in nearly all the applications of the literature) may inherently be written as a linear combination of those structural displacements where a load is applied (and in fact the loads themselves are the coefficients of such linear combination as shown in the sequel of the paper). Any optimal design pursued herein is the one that minimizes a proper measure of "amplitude" of the transfer matrix G . From an engineering perspective, minimizing such amplitude amounts to minimize the effects of the loads onto the structure. It is therefore clear that a proper "tool" to assess the amplitude of the transfer matrix is the key ingredient to fully define the proposed approach. The Singular Value Decomposition (SVD) (Brunton and Kutz 2019;Strang 2019) is shown to be such a tool and more in that not only it does allow to compute a few matrix norms that represent the needed amplitude but also shed some light onto the system response in static and dynamic regime as well. The fact that such norms are uniquely defined in terms of singular values of G simplifies dramatically the computation of the sensitivities that are reduced to the evaluation of the gradients of the singular values themselves with respect to the design variables.
The paper is organized as follows. Section 2 introduces the properties of the singular value decomposition focusing on its engineering interpretation. The fact that input and output singular vectors are mapped onto each other with a gain that is the associated singular value is the reason why the proposed approach is in a sense the most natural to address static and dynamic topology optimization problems, including the celebrated minimum compliance one. The section is concluded by an explicit derivation of the transfer matrix in the frequent case of a displacement-based FEM discretization. The topology optimization problem is presented in Sect. 3 along with a semi-analytic derivation of the sensitivities of the transfer matrix entries and of the singular values with respect to the design parameters. Such sensitivities are needed to apply any gradient-based minimization scheme including the method of moving asymptotes (MMA) (Svanberg 1987) that is used herein in a Matlab version made available by Krister Svanberg. The specific choice of the output vector that fully defines the objectives of the design is made in Sect. 4 wherein attention is paid to the compliance or better a few generalized versions of it. The extension to the dynamic problem is analyzed in Sect. 5. The focus therein is not on the (minimal) modifications needed to handle the dynamic case but is conversely on the unifying aspects that static and dynamic problems share. A few convenient computational remedies that take care of the lengthy processes behind the expensive evaluation of system matrix norms are highlighted and framed within the Matlab environment. Section 6 presents the results of a few numerical investigations along with comparisons with similar designs obtained via more classical methods (when available). The focus is on reporting general lessons that may be learned by proper comparisons between all the proposed norms in terms of effectiveness and robustness. Concluding remarks and need for ongoing investigations are left to Sect. 7, whereas the full Matlab code that solves the static topology optimization problem is given in online Appendix 1. The code is provided for completeness and is open for extensions and improvements.

General properties of the SVD
• U and V are unitary l × l and m × m matrices, respectively; • is a l × m matrix that admits the representations where in which r is the rank of G and A few developments below are based on the fact that the unitary matrices U and V form orthonormal bases for the output and input spaces of G , respectively.

The main property
Equation (2) provides the definition itself of singular value decomposition but should be further elaborated to gain a deeper understanding of the proposed design approach. To this goal, both sides of Eq.
(2) are right-multiplied by V so as to obtain (thanks to the unitary property of V) that componentwise reads Equation (5) is the first of the two milestones on which the proposed design approach rests (the second one is concerned with the concept of matrix norms based on singular values that is tackled in the next paragraph). The column vectors of U and V , say u k and v k , respectively, represent the output and input directions of the system and are dual to each other, i.e., if an input vector (load) in the direction of v k is considered, the output shall be in the direction of the associated singular vector u k . Since ||v k || 2 = ||u k || 2 = 1 , one has From an engineering point of view, Eq. (6) is of the utmost importance for two reasons: -The kth singular value k is the gain of the matrix G in the direction of the k− th input-output singular pair; -Being singular values sorted in descending order, see Eq. (4), the triple ( 1 , v 1 , u 1 ), respectively, define the maximum gain of the response matrix G and the associated input and output directions.

Singular values and matrix norms
Paragraph 2.1 has introduced the singular value decomposition of the response matrix as an engineering tool for a deeper understanding of the input-output mapping of the structure to be designed. Furthermore, singular values were shown to be the structural gains coupling each singular vector pair. Given this scenario, a natural goal function, say F , should enjoy the following two properties: 1. F should be a proper norm of the response matrix G (and such norm should preferably be induced by an associated vector norm); 2. The norm F should be uniquely defined by the singular values (gains) k .
Three possible such norms are as follows: • ||G|| 2 = 1 (the so-called 2-norm) that is in fact the supremum of the ratio in (6). The 2-norm satisfies the multiplicative property which is a desirable property that shows that the 2-norm is a vector-induced one. As shown later in much detail, it is appropriate to note that in the dynamic case singular values depend on the frequency , i.e., 1 = 1 ( ) , and the induced norm is defined as sup = 1 ( ) . Since a sup operator is involved, in dynamics such norm shall be named H ∞ norm but one should always recall that the underlying norm is a 2-norm; • ||G|| N = 1 + 2 + … + r that is often referred to as the Nuclear norm. All plant's directions concur to the definition of the goal function and therefore the frequency range wherein the design is effective is expected to be wider than the H ∞ -norm case that however performs better as far as maximum peak reduction is concerned; • ||G|| F = √ 2 1 + 2 2 + … 2 r , i.e., the so-called Frobenius norm. Although neither is the Frobenius norm an induced one nor it enjoys the multiplicative property (7), numerical results to follow show that the Frobenius norm works properly and is a viable choice when the primary goal is the robustness of the response.
To conclude this section with, the following key features of the singular value decomposition are worth being emphasized: -the singular values give direct info on the input/output gains of the system; Page 4 of 27 -the structure directions U k and V k are mutually orthogonal and represent the input-output directions of the system; -unlike eigenvalue-based approaches, the SVD applies inherently to non-square structures. Most response matrices G are in fact non-square since the number of inputs and outputs are usually different from each other. (There is actually no reason why they should be the same).

The steady-state input/output response matrix G
Aim of this paragraph is to derive the abstract framework for the static state and output equations that fits the proposed SVD-based design approach. To start off, equations governing the static response of the structure under optimization are written. For the sake of simplicity, reference is made hereinafter to a displacement-based finite element discretization, but mixed stress-displacement approaches based on the Hellinger-Reissner functional or three-field (stress, strain, displacement) formulations arising from the Hu-Washizu functional might be considered as well. A classical displacement-based equilibrium equation reads where K is the n × n stiffness matrix, x is the n × 1 column vector of the unknown displacements, B is an n × r load matrix, and f is a r × 1 Boolean column vector that distributes the loads to the appropriate degrees of freedom. It should be noticed that constrained displacements and dual reactive forces have been ruled out from Eq. (8) so that n denotes the number of free degrees of freedom of the problem. The output equation is then used to define the goals of the design. Each row of the n o × 1 output vector z represents a physical quantity of interest for the designer. When compliance minimization is pursued n o = 1 so that z is actually a scalar quantity. The output equation is eventually written as follows: where C is the n o × n output matrix that is responsible for the definition of the quantities of interest for the designer.
Upon combining Eqs. (8) and (9), i.e., one gets the input-output mapping that reads that is given a classic system representation in Fig. 1 where the dependence of G on the design variable vector p has been made explicit.
The n o × r steady-state (static) transfer matrix eventually reads where one should notice that in the frequent case of a single load case ( r = 1 ) and a single scalar output ( n o = 1 ), as it happens, e.g., for compliance minimization of a structure acted upon by a single load, G is actually a scalar quantity.

G-based topology optimization
The typical G-based optimal design problem eventually reads where -(12) 1 defines the function F(p) to be minimized. The matrix norm being used for each specific problem should be specified (see paragraph 2.2); -(12) 2 is the definition itself of the input/output transfer matrix wherein the dependence on the vector of the design variables p has been made explicit; -(12) 3 is the standard constraint on the maximum volume of the structure. From a computational viewpoint, it may be useful to impose also a lower bound to the overall volume, i.e., Although not compulsory, in the very first iterations of the optimization process a lower bound on the volume

Fig. 1
Block representation of the input-output mapping may help the numerical solver avoid unfeasible "superlight" solutions in the design space; -The inequalities on the design variables in (12) 4 are to be intended componentwise. As usual, each design variable is bounded between zero (void) and one (full material) within a standard SIMP idealization (Bendsøe and Sigmund 1999).
The optimization problem (12) is solved numerically using the method of moving asymptotes (MMA) (Svanberg 1987). To this goal the user must provide the current values of the goal function and the constraints, along with their sensitivities with respect to the design variables. When the objective function to be minimized is a matrix norm that depends on its singular values, using the chain rule one should therefore provide a semi-analytical computation of the sensitivity of the transfer matrix and of the singular values with respect to any design variable, say p k . This is done below and implemented in the Matlab code in online Appendix 1.

Sensitivity of the input/output transfer matrix
In the most general case assume the input/output matrix G depends on the design variable vector p through all its three factors, i.e., In fact, when, for example, seismic forces are considered, one may show that the matrix B does depend on the design variable vector p (Venini and Ceresa 2018) and that is why the dependence of all the factors in Eq. (13) on p should in general be considered. Upon recalling that the derivative of the inverse reads using the product rule, the derivative of G(p) with respect to the kth component p k of the design variable vector writes where the derivative of the stiffness matrix K(p) with respect to the design variable p k that appears in Eq. (14) is computed semianalytically (Andreassen et al. 2011). When the dependence of G on p is through K only, Eq. (15) reduces to

Sensitivity of the singular values
Under certain hypotheses of regularity (that we certainly supposed to be met hereinafter) direct differentiation of Eq.
(2) allows to arrive at the following formula for the sensitivity of the generic h-th singular value h (see also Strang 2019): where R indicates the real part of a (possibly) complex number.

Sensitivity of the norm of the input/output transfer matrix
Since 2, Nuclear and Frobenius norms of G are uniquely defined by its singular values, Eq. (17) coupled to a proper usage of the chain rule allows the computation of any norm sensitivity in (12) 1 . The sensitivity of any of the three norms of interest eventually reads where the second factor at the right-hand-side is given in Eq.
(17), whereas the first factor is as follows:

Compliance definition (revisited and extended) and more
As already observed most of the times topology optimization aims at the minimization of the compliance. Therefore, even though the proposed approach may handle many different goal functions, a specific derivation of the (16) compliance minimization problem seems to be appropriate. For the sake of simplicity, a lumped parameters system is referred to next. However, the adopted formulation is general in that it applies to any system such that x is the state vector grouping all the generalized displacements emerging from a suitable finite element discretization; -s ⊂ x is the subset of the displacements that is loaded by the force vector f that is the dual of s in the virtual work sense.

Notations
Reference is made to Fig. 2 where the state vector reads The optimization problem associated to the system in Fig. 2 is far from being unique in that different choices of load cases and combinations on the input side, see Eq. (8), and design objectives on the output side give rise to a wide variety of engineering problems of interest as discussed below. Each problem fits the framework of Eq. (12) and is fully characterized by the stiffness matrix K , the load matrix B , and the output matrix C . For the case at hand n = 3 and the 3 × 3 stiffness matrix is shared by all possible problems. The designer should conversely select the two matrices B and C to derive the optimization setting that fits better the actual engineering problem to be solved.

The SISO case: single load case and standard compliance minimization
In this case, i.e., the standard compliance minimization in the presence of a single load case, n o = 1 , r = 1 and one writes where it should be noticed that the output vector z reduces to a scalar and reads: i.e., the classical compliance (Andreassen et al. 2011). In this SISO case G is a scalar as is the "matrix" of singular values that reads = 1 . It should be noted that in the SISO case the three cited norms (2, Nuclear and Frobenius) coincide.

The MISO case: two load cases and standard compliance minimization
This case is non standard in that r = 2 , i.e., two independent load cases are considered, whereas n o = 1 , i.e., the output is a scalar compliance. One may write where the output vector is formally the same as the SISO one but now the overall compliance C is the sum of two independent ones, say C 1 and C 2 , respectively, due to the loads f 1 and f 2 , i.e., This resembles the multiple load formulation of Sect.5.2 in the celebrated paper (Andreassen et al. 2011) by Sigmund and co-workers. A further property enjoyed by the MISO case should be highlighted. Having two inputs and a single output, the singular input matrix V and the singular value matrix may be, respectively, written as The only non-vanishing singular value 1 gives the plant gain, whereas the associated first right (input) singular vector v 1 gives the factors of the linear combination of the two loads f 1 and f 2 that stresses the structure the most, i.e., On the opposite side, the second right singular vector v 2 is associated to the vanishing singular value 2 = 0 . Therefore, to which a null compliance is associated. The crucial role played by (26) as far as the optimal topology is concerned should be highlighted: as long as the two loads f 1 and f 2 are independent of one another, i.e., when they belong to different load cases so that r = 2 , the proposed approach provides the designer with an optimal solution that is based on the automatic selection of the worst possible load combination that is given by the components of v 1 and characterized by a gain that is 1 . No role is conversely played by the load combination (27). Being the system a single output one, i.e., (25), the rank of G is r = 1 and therefore the plant direction v 2 is annihilated or, in engineering words, the relevant compliance is zero whenever the two loads f 1 and f 2 are combined according to (27).

The SIMO case: single load case and minimization of the compliance vector norm
The peculiarity of this (and next) case is its output z . The compliance associated to each loaded point, say C k = f k u k , is no longer added to define a unique scalar compliance C = ∑ k C k , but it conversely defines an independent entry of the compliance vector z . For the problem at hand r = 1 and n o = 2 and one writes such that the output vector reads The rectangular matrix of singular values reads that means that C 1 and C 2 are automatically combined according to the components of the first output vector u 1 , i.e., and such combination is the one that determines the final optimal topology to which the maximum gain 1 is associated. On the other side, the combination is associated to a null gain 2 = 0 and is therefore annihilated by the computations. As a matter of fact (31) is the SIMO dual of the MISO input combination (26) in that they both are the combinations, input and output, respectively, that govern and determine the optimal topologies. Likewise, (32) and (27) are dual to each other in that none of them plays any role toward the computation of the final optimal topologies being associated to 2 = 0.

The MIMO general case: two load cases and minimization of the compliance vector norm
This is the way more general formulation. One has two inputs and two outputs, i.e., r = 2 and n o = 2 , thus generating a 2 × 2 transfer matrix G . As for load and output matrices one may write the output vector z being the same as in Eq. (29). In this case ( 2 = r > 1 and 2 = n o > 1 ), two non-vanishing singular values show up and one may write where 1 and 2 are, respectively, associated to the structure directions u 1 and u 2 that maximize and minimize the plant gains. In the MIMO case, a comparison between the choice of the 2-norm on the one side versus the Nuclear and Frobenius norms on the other is appropriate. Being r = n o = 2 > 1 , no plant direction is annihilated (unlike previous cases) and both the triples ( 1 , u 1 , v 1 ) and ( 2 , u 2 , v 2 ) contribute to the output response. Being however 1 > 2 , the first triple plays a more important role than the second. That said, the 2-norm is focused on the first triple only, i.e., the one that maximizes the plant gain, whereas using the nuclear and Frobenius norms all plant directions enter the formulation, including those characterized by smaller gains ( 2 for the case at hand).

Dynamic topology optimization fits the proposed formulation
By definition itself, the proposed topology optimization approach of Eq. (12) may be applied to any system whenever the following properties are met: When these two properties hold true, the proposed formulation is abstract enough to allow its applicability, no matter the specific problem at hand and in particular no matter whether the system is static or dynamic. For the latter to fit the proposed formulation, a Laplace (Fourier) domain setting is adopted wherein the transfer matrix G is readily available in an overall algebraic framework. In the dynamic case, though, the transfer matrix is no longer constant but depends on the Laplace variable s and eventually on the frequency (and this matter shall affect the definitions themselves of matrix norms).

The dynamic transfer function
The more convenient way to arrive at the Laplace-domain formulation of the dynamic problem is to start with the time-domain state space form and take its Laplace transform. The state space form reads where a superposed dot indicates time differentiation and in which -K , D , M , and I are the stiffness, damping, mass, and identity matrices, respectively. In what follows a Rayleigh damping is considered such that where and are positive constants (small enough not to overdamp the system). However, the formulation is not limited to Rayleigh damping and any viscous damping could be adopted; -x is the state vector that piles the displacement and velocity vectors u and u , respectively; -z is the output vector that is defined through a proper choice of the output matrix C.
Taking the Laplace transform of (35) 1 and (35) 2 and ruling out the Laplace transform of the state vector x one gets the Laplace-domain input/output relationship and the desired transfer matrix G(s) , i.e., Upon choosing s = i one gets the frequency response matrix that is the dynamic counterpart of the static one in (12) 2 : The dependence of G(i ;p) on the frequency ∈ (0, ∞) calls for a generalization of the three norms previously introduced as discussed below.

Matrix norm definitions in the dynamic case
The norms in Sect. 2.2 are extended to the dynamic case in which G is no longer constant but depends on the frequency . This is done through an operation that is conceptually trivial but computationally extremely expensive: the supremum over the frequency axis in the H ∞ and Nuclear norm cases, a suitable integral average in the H 2 case. One may write -H ∞ norm: ||G|| ∞ = sup 1 ( ) where again the slight change of notation should be noticed. In statics this was named a 2-norm in view of the fact that it was defined as a ratio between the 2-norm of the output and input vectors.
In dynamics, such norm is named H ∞ since a supremum is taken over the frequency axis; -nuclear norm: ||G|| N = sup 1 ( ) + … + r ( ); -the H 2 -norm is the generalization of the algebraic Frobenius norm to dynamical systems. Unlike the H ∞ -norm, however, such generalization is not performed using a sup operator, such as but via an integral average that reads However, Eq. (40) is neither practical for computing the H 2 norm nor its sensitivity. A more viable approach is then recalled in Sect. 5.3 below. (38) Page 9 of 27 85

Computing the H 2 norm and its sensitivity
The first step to compute ||G|| 2 in the dynamic case is to evaluate the contrallability Gramian W that solves the generalized Lyapunov equation: Then, one may show that the following relation holds true (Francis 1987):

Semi-analytic H 2 -norm sensitivity
As for the derivative of ||G|| 2 with respect to the design parameter p k , one may write by direct differentiation of (42): in which the derivative of the Gramian W with respect to p k may be computed by solving the generalized Lyapunov equation one gets by direct differentiation of (41): where a superposed prime denotes differentiation with respect to p k and the pseudo-load S is given as follows:

Finite difference H 2 -norm sensitivity
The semi-analytic sensitivity of Eq. (43) calls for the solution of an additional generalized Lyapunov equation, the pseudo-load of which is quite a complex one, i.e., S in (45). A fast solver for the generalized Lyapunov equation was coded in the Matlab add-on Saak et al. (2022) that performs at its best when the B matrix is a "thin" rectangular one, a property that is not enjoyed by S in (45). It is therefore more convenient to compute the H 2 -norm sensitivity using a finite difference approximation that also involves the solution of a generalized Lyapunov equation. The B matrix of the latter, though, is the same as the one of the actual system and is therefore ideally suited for the code in Saak et al. (2022) to work full throttle.

The dynamic topology optimization problem
The dynamic topology optimization problem eventually writes that is formally the same as the static one in (12), provided the respective definitions of matrix norms (Sect 2.2 versus Sect 5.2) and transfer matrices (Eq. 11 versus Eq. 39). A few considerations on how to compute the norms and their sensitivities when the frequency response matrix depends on are in order. The scope here is to show how these norms can be computed using, for example, Matlab, whereas details concerning the algorithms themselves are beyond the scopes of this paper and may be found in Bruisma and Steinbuch (1990) among others.

On the numerical computation of ||G(i!, p)||
Computing any of the three norms of G(i ;p) (for fixed p ) represents a formidable computational task. The H ∞ norm is by far the one that has received more attention by the scientific community, mainly because of its relevance in the framework of robust control of closed-loop systems (Francis 1987). The Matlab function norm (or better its version made available in the add-on sss Toolbox (MORLab 2017) that works on sparse matrices) may be used for this purpose: provided as inputs a suitable version of the transfer function G and the flag , the function norm returns the H ∞ -norm of G and the frequency at which the norm is attained. The latter frequency is needed to compute the sensitivities as shown below. Using the same function with the flag 2, the Frobenius norm is retrieved. However, when the system size is 10 5 or more as is the case of the examples below, all available norm functions become dramatically slow and the very many calls required to achieve convergence within an optimization session make their usage not convenient if not impractical. The way out adopted for the simulations to come is a sub-optimal strategy described below that applies to the Nuclear norm as well for which ad hoc methods are not available.

A viable sub-optimal approach for computing the H ∞ and the Nuclear norms of a large scale dynamic system
Given a dynamic system with transfer matrix G(i ) , the proposed method for a sub-optimal computation of ||G(i )|| ∞ and ||G(i )|| N works as follows: 1. the frequency range of interest, say I = ( 1 , 2 ) , is determined and sampled (uniformly or not); 2. a call to the Matlab function sigma is made that returns all the singular values of the dynamic system for all frequencies in I ; 3. ||G(i )|| ∞ and ||G(i )|| N are, respectively, computed as follows: The method is only sub-optimal in that the first singular value and the sum of all of them might assume larger values at frequency points that do not belong to the set I , especially if the adopted partition of the frequency axis is a coarse one. One may of course refine the approach using an adaptive sampling strategy that would allow a closer and closer computation of the actual value of the norm. For the sake of this paper, this option is not however exploited but is left for near-future extensions.

On the computation of the sensitivity of ||G(i!, p)||
If one compares the static and dynamic definitions of G that are, respectively, given in Eqs. (12) 2 and (46) 2 , it is straightforward to realize that the static sensitivity formulas in (15) and (16) still hold true in the dynamic case provided that -K is substituted by i E − A; -in such formulas, is no longer a variable but takes on the value at which the norm is attained.

Numerical results
A few representative numerical results are presented hereafter that are cast in the framework of Sects. 3 and 5 for the static and dynamic cases, respectively.

Software and third party add-ons
The codes with which all the numerical results were obtained are written in Matlab. The full code that solves the static problem is presented in online Appendix 1 and requires the Matlab implementation of the Method of Moving Asymptotes (MMA) made available by its creator Krister Svanberg, Svanberg (1987). The following add-ons were used for the dynamic code: is an ideal tool in this respect, e.g., for solving the generalized matrix Lyapunov Equation (41).

Preliminary considerations
The four possible formulations (SISO, MISO, SIMO, and MIMO) are exploited along with the three proposed norms (2, Nuclear and Frobenius). It should be outlined that in the first three cases (SISO, MISO, and SIMO) r = 1 , i.e., a single singular value is non-trivial so that the three norms happen to be the same and so do the relevant optimal topologies. Therefore, for each of the cases SISO, MISO, and SIMO a single optimal topology shall be shown that is generated using any of the three norms. In the MIMO case, though, two singular values are non-trivial and therefore the optimal designs do depend on the norm being used. Therefore, for the latter MIMO case each design is presented.

General data
The geometry and the loads of the static topology optimization problem are shown in Fig. 3. Thanks to the symmetry of the structure the half-model in Fig. 3 is used in the computations using the following numerical data: -Young modulus: E 0 = 1 , E min = 10 −9 , -Poisson ratio: = 0.3, -Number of elements in vertical (Y)-direction: 2 6 ; -Number of elements in horizontal (X)-direction: 3 × 2 6 ; -Side length: L = 2 6 (by this choice a uniform mesh of square elements of unit size is generated); -P = Q = 1 (the load time dependence in Fig. 3 is for use in the dynamic example below); -Volume fraction = 0.5; -Filter type: density, see Andreassen et al. (2011); -SIMP penalization factor: 3.

SISO formulation
Using a SISO formulation the input is given by the two forces P and Q that belong to the same load case and the output is the overall compliance C . As shown in Table 1, the norms at convergence coincide (as they should), whereas the computational times do not. The CPU time is limited and smaller than two minutes in all cases, see the convergence curves in Fig. 4. The optimal topology is shown in Fig. 5 that is the same as the one given by the 88-line code by Sigmund and co-workers Andreassen et al. (2011). This does not come as a surprise since the problem being solved is the same, even though the underlying formulation is different. As for the input and output plant directions, in the case of a SISO representation both V and U degenerate to the (scalar) unity, whereas the singular value is the norm itself.

MISO formulation
In the MISO formulation, the two loads P and Q are allowed to vary independently of each other within two different load cases. However, a single output compliance is considered that is the sum of the two compliances. Table 2 shows norm values and CPU times needed to perform the computations. The convergence curves are shown in Fig. 10 that exhibits a "nearly" monotonic convergence for all the three norms. The optimal topology (shared by all the three norms) is shown in Fig. 7. A look at the topologies in Figs. 5 and 7 shows that, though similar, the two designs are different. One may wonder whether this represents a violation of the superposition principle according to which as long as the response is linear the output of the sum ( C due to P and Q acting upon the structure simultaneously) is the sum of the outputs ( C = C P + C Q ). Of course, no violations, since the structure is not given but being defined iteratively, the updates on the design variables at each iteration need not to be the same and so are the final designs (Fig. 6) .
The input plant direction is given by the singular input vector at convergence that reads However, for MISO systems, the concept of plant direction is not as important as one may show that the gain neither depends on the input magnitude nor on its direction (Skogestad and Postlethwaite 1996), i.e., V . As shown below, the concept of input direction becomes conversely crucial in the case of MIMO systems.

SIMO formulation
In the SIMO formulation, the two forces P and Q belong to the same load case but the overall compliance is split to form a 2D output vector z = C 1 C 2 T . Table 3 shows norm values and CPU time needed to achieve convergence. Figures 8 and 9 show the convergence paths for the three norms and the optimal topology, respectively.
The optimal topology happens to be the same as the MISO case and interestingly enough the singular output vector reads that is the same as the input vector of the SIMO case.

MIMO formulation
As anticipated, the MIMO case is the more general one in that not only are the two forces P and Q free to vary independently of each other within two load cases, but also each compliance C P and C Q represents an entry of the 2D output vector z . The concepts of plant directions, orthogonality of input and output singular vectors, and singular values seen as plant gains assume the deepest meaning in the MIMO case. Norm values and CPU time needed to achieve convergence are displayed in Table 4.
The optimal topologies obtained using the three norms are shown in Fig. 11. For the sake of completeness, Fig. 12 shows the optimal topology one gets using Sigmund's 88-line code with the multiload option. The four designs are similar but different being based on different optimality conditions. The triples (V, U, ) are reported next for the four designs at convergence:     expected because of the cantilevered nature of the beam and the position of the two loads (P and Q, respectively, inside and outside the beam span). As for the performance of the proposed designs, an insight may be gathered by looking at Table 5 that grades the presented designs in terms of the numerical values of the adopted norms. The following considerations are worth being made: -as expected each design wins its own competition, i.e., any norm is minimized by the formulation that explicitly aimed at its minimization (the relevant values are written in bold on the diagonal of Table 5); -minimizing the 2-norm, i.e., the maximum gain 1 of the system in the direction of the first input and output singular vectors, causes a high penalization on the other plant directions. One may in fact see that for the three cases Nuclear, Frobenius, and Sigmund, the second (and last) singular value is, respectively, equal to 38.22, 41.57, and 38.06, i.e., values that are smaller than half the first gain 1 . Not so for the 2-norm case for which 2 = 54.72 , i.e., a remarkably larger value. As a general guideline, the lesson learned is that using a norm that involves the maximum singular value only is computationally more efficient and inherently leads to the reduction of the maximum gain (which is an amazing result). However, neglecting the effect of smaller singular values may lead to a system response that is poorer overall, especially when 2 is close to 1 ; -for each of the four designs, a measure of overall performance (the smaller the better) may be computed as the SRSS (square root of the sum of the squares) of the differences between scored norm values and best norm values (eventually scored by other approaches). One gets the results shown in Table 6 that again shows that the 2-norm is an extreme one: the maximum gain is reduced the most among all possible design choices but the price    to be paid is a weaker behavior with respect to other possible merit functions. Eventually, Sigmund's code design (that is not based on an explicit norm minimization as other designs do) is shown to be quite a robust one in that it performs nicely with respect to all the proposed norms.

Preliminary considerations
The five-story high-rise building shown in Fig. 13 is investigated next. The objective is to arrive at an optimal bracing structure capable to withstand horizontal environmental actions. This kind of application of topology optimization has been investigated for some fifteen years (see Beghini et al. 2012 as one of the most successful strategies within this framework). This problem is properly suited to test the proposed approach in that its multi-output version is a five outputs challenging one. The following four different designs are pursued: 1 Convergence curves for all the designs are displayed in Fig. 14 that show a steep monotonic behavior of the MMA-based approach. All the optimal topologies are shown in Fig. 15. The SISO-based bracing system is shown to be quite a stiff one at the base (and is so up to roughly mid-rise). MIMO designs are in general characterized by a more homogeneous material distribution as far as the thickness of the diagonal braced is concerned. From a manufacturability standpoint, the MIMO design using a nuclear norm-based goal function appears to be superior with respect to other designs in that a lower and, so to speak, "cleaner" bracing layout shows up. The performance of three MIMO designs is assessed in Table 7 where all the three norms are computed.

Preliminary considerations
Attention is focused on the MIMO formulation only that is the one that exhibits a more general and structurally richer behavior. As for the norms, the H ∞ and the Nuclear norms are considered since they depend on local evaluations of the singular values in such a way that the sub-optimal approach introduced in Sect. 5.4.2 may be used. On the other side, though conceptually applicable, the usage of the H 2 -norm has not been exploited because of the prohibitive computational cost of the relevant sensitivities, no matter the function or approximation being used.

General data
Geometry and loads are the same as Problem 1. As for the Young modulus E 0 = 1000 and E min = 10 −9 are used, whereas the mass density parameters are chosen as 0 = 10 −3 and min = 10 −6 and a Rayleigh damping matrix is introduced that reads The number of elements in vertical and horizontal directions are given as 2 5 and 3 × 2 5 , respectively.

H ∞ -norm formulation
The optimal topology one gets by minimizing the H ∞ -norm of the MIMO transfer matrix G(i ) is shown in Fig. 16, whereas Fig. 17 shows the path of the MMA scheme toward convergence. The chattering pattern displayed by the convergence curve is not due to numerical difficulties encountered along the path by the optimization routine. It is conversely due to the fact that, as explained below in more detail, the first singular value 1 ( ) exhibits two very close points that are competing to realize the overall maximum, i.e., the H ∞ -norm. The design switches back and forth between the two points along an overall converging path. Eventually the two designs coalesce and convergence is achieved. Figures 18 and 19 shade some light on the complexity of the dynamic response of a Multi-Input Multi-Output system (even the simplest of such type, i.e., a 2-input 2-output structure as the one under investigation) and allow a deeper and rationale understanding of the overall optimization process. Remarkable features emerging from Fig. 18 are as follows: -Toward the evaluation of ||G(i )|| ∞ , i.e., sup 1 ( ) , a competition between two frequencies takes place to assess the one at which the norm itself is attained. This should be considered a general outcome of minimizing the H ∞ -norm. In fact such an approach is focused on minimizing the supremum of 1 (the maximum of the first singular value) over the frequency axis. Design changes are therefore sharply localized in the frequency range where the current maximum is being attained, whereas no care is taken of the dynamic response at frequencies far apart at which the magnitude of the response is allowed to increase. When such an increase leads to a response amplitude that falls closer and closer to the one where the H ∞ -norm is currently attained the competition starts. In general, one can conclude that designs based on the minimization of the H ∞ -norm are likely to be multimodal ones as a result of this inherent feature; (53) D = 0.04M + 0.03K. -it is well known that classical eigenvalue optimization may be affected by the eigenvalue crossing phenomenon, i.e., a crossing between different eigenvalue trajectories in the design space that leads to the annoying non-differentiability of the min − max eigenvalue problem. A similar behavior shows up in Fig. 18 where a singular value crossing takes place at 4.0 rad/s. Such crossing takes place at a minimum point so that the proposed procedure is not affected by the phenomenon.
However this may not always be the case and an issue may arise in general so that further investigations in the near future are planned on this topic.    Switching to the amplitude of the frequency response, reference is made to Fig. 19 where the points in Table 8 are highlighted. First of all one should observe that the two points P 12 and P 21 govern the mutual compliance, i.e., the (amplitude of the) displacement of B when the loaded point is A and viceversa. The fact that P 12 = P 21 is a version of the Maxwell-Betti reciprocal work theorem applied to the dynamic case in a harmonic load framework. Interestingly enough, not only are the two amplitudes the same, but also the frequencies at which maximum points are attained coincide. As for the overall system dynamics, the singular value crossing at = 4 rad/s and the competition between the two maximum points find their motivation and confirmation in Fig. 19 where one may see that the amplitude of the compliance of the two loads (points P 11 and P 22 ) is quite similar as are the relevant frequencies P 11 = 3.76 rad/s and P 22 = 4.49 rad/s. From a quantitative point of view, the singular value decomposition of G(i ) for all the frequencies above read as follows: 1.

Nuclear norm formulation
The optimal topology one gets by minimizing the Nuclear norm of the MIMO transfer function G(i ) and the    Fig. 21 no oscillations show up and this is a first indicator of a structural response dominated by a single mode. This is confirmed by the singular value curves displayed in Fig. 22 wherein an innocuous low-amplitude crossing takes place at about = 5 rad/s. In fact the second gain (maximum of the second singular value) is smaller that half the larger gain that takes place at = 3.98 rad/s that therefore governs the entire system dynamics. Figure 22 shows the variations of the singular values 1 ( ) and 2 ( ) and their sum 1 ( ) + 2 ( ) with respect to the frequency .
Outstanding features emerging from Fig. 22 are as follows: -the Nuclear norm, i.e., the maximum value of 1 ( ) + 2 ( ) is 0.8789 and is attained at = 3.98 . Using a goal function that depends on both the singular values causes design updates to spread their effect over a wider frequency range, at least as wide as ( P 11 , P 22 ) , see also Table 9. This seems to be the main reason why the amplitude of the second singular value is kept small and a single-mode design is ended up with, unlike the H ∞ -norm-based design that was a multi-modal one; -Not only is the dynamic response dominated by a single mode but also a similar behavior is exhibited by the two inputs. Upon comparing the upper and lower amplitude curves in Fig. 23, one concludes that the plant is dominated by the first input, i.e., P, whereas the second input Q has a nearly negligible effect.
To summarize with, (at least for the problem at hand) using the nuclear norm leads to a single-mode design wherein one load determines much of the structural response, whereas the other is nearly annihilated by the design itself. These features are not shared by the H ∞ -norm-based design that puts all its efforts against the overall peak response.
The single modality of the structural response is further confirmed by the singular value decomposition of G for = P 11 , = P 22 , and = P 12 . In all cases, the second (and last) gain 2 is much smaller than the largest gain 1 as shown below:

Conclusions and need for further investigations
A novel unifying approach has been presented that allows to handle static and dynamic topology optimization using the very same formulation (and computer code). The singular value decomposition of the input-output transfer matrix G along with the selection of suitable matrix norms (defined uniquely in terms of its singular values) to be used as goal functions represent the key ingredients of the proposed approach wherein static and dynamic frameworks are identical. Computing G , its norms and its sensitivities in the static and dynamic cases are the main focus of the paper that is completed by in-depth numerical investigations wherein general features of the resulting designs are highlighted. The full Matlab code that solves the static problem is provided in online Appendix 1 and may be readily extended to cover the dynamic case by properly implementing the details given in Sect. 5. Ongoing investigations are concerned with the following topics that are mainly focused on the improvement of the dynamic topology optimization approach: • The number of states of the dynamic systems investigated herein is about twenty thousands, whereas just few poles govern the dynamic response of the system. Using reduced order models seems therefore the more natural way to cope with dynamic accuracy and reduction of CPU time; • Using reduced order models opens the way to the potential issue of computing the sensitivities. In fact reduced order models are in general black boxes that do not allow design variables at the element level to hold their physical meaning. A possible way out is the adoption of parametric reduced order models that are currently being investigated; • The H ∞ -norm approach was shown to lead to multimodal solutions characterized by a competition between candidate design points at different frequencies. The numerical minimization scheme is therefore likely to be trapped in between the two points exhibiting an overall chattering behavior. The adoption of proper filters that allow a fast choice among the two points is the remedy being exploited to solve such problem; (70) 12 = 0.641 0 0 0.088 , (71) V 12 = −0.845 + 0.000i − 0.535 + 0.000i 0.534 − 0.027i − 0.844 + 0.043i .

Fig. 23
Nuclear norm formulation -Frequency response amplitudes-Top: from P (input 1) to outputs 1 ( u A ) and 2 ( u B )-Bottom: from Q (input 2) to outputs 1 ( u A ) and 2 ( u B ) • There are such applications as, for example, building engineering wherein structures should be designed to withstand static and dynamic load combinations at the same time. The former are often determined by permanent and variable vertical loads, whereas earthquakes represent the driving loads of the latter. It is therefore planned to extend the proposed approach using goal functions that are convex combinations of the purely static and dynamic ones; • Finally a few words of the robustness issue (both in statics and dynamics). It is planned to adopt some nonprobabilistic uncertainty models similar to what is done in the framework of robust control of uncertain systems. By doing so, one operates on a class of systems (and not on a single one) in a worst-case scenario framework; To conclude with, it is planned to use some Application Programming Interface (API) to link the Matlab code to a general purpose finite element software so as to apply the proposed approach to (nearly) real-world structural and industrial applications.