Integrated topology and controller optimization using the Nyquist curve

The design of high-performance mechatronic systems is very challenging, as it requires delicate balancing of system dynamics, the controller, and their closed-loop interaction. Topology optimization provides an automated way to obtain systems with superior performance, although extension to simultaneous optimization of both topology and controller has been limited. To allow for topology optimization of mechatronic systems for closed-loop performance, stability, and disturbance rejection (i.e. modulus margin), we introduce local approximations of the Nyquist curve using circles. These circular approximations enable simple geometrical constraints on the shape of the Nyquist curve, which is used to characterize the closed-loop performance. Additionally, a computationally efficient robust formulation is proposed for topology optimization of dynamic systems. Based on approximation of eigenmodes for perturbed designs, their dynamics can be described with sufficient accuracy for optimization, while preventing the usual threefold increase of additional computational effort. The designs optimized using the integrated approach have significantly better performance (up to 350% in terms of bandwidth) than sequentially optimized systems, where eigenfrequencies are first maximized and then the controller is tuned. The proposed approach enables new directions of integrated (topology) optimization, with effective control over the Nyquist curve and efficient implementation of the robust formulation.


Integrated controller-structure optimization
Many high-tech applications require positioning at both high accuracy and high speed, for which motion systems are used. These are, for instance, used in semiconductor equipment, microscopy, robotics, and medical devices (Munnig Schmidt et al. 2011;Oomen 2018). The required speed and accuracy in these positioning systems is achieved by feedback control. In the quest for more extreme performance, the design of motion systems poses a significant challenge.
The final performance and accuracy of such systems heavily depend on system dynamics, the controller, and the (closed-loop) interaction between the two (i.e. mechatronics). Various complex design problems have been effectively addressed by topology optimization in recent years, and the need exists to also apply it to the design of motion systems. Although optimization is frequently used in the design of feedback controlled systems, it is mostly applied in a sequential manner. First, the structure is designed, e.g. for maximum eigenfrequencies using topology optimization (Ma et al. 1995;Delissen et al. 2022), after which a controller is tuned for this structure that achieves system requirements, such as high bandwidth, closed-loop stability, and disturbance rejection (Munnig Schmidt et al. 2011). However, this approach usually leads to sub-optimal system performance. High eigenfrequencies are often a characteristic of good system performance, but it does not mean that higher eigenfrequencies always result in a higher bandwidth. Therefore, for superior performance of the combined system, an integrated approach is required (Fathy et al. 2001; van der 1 3 80 Page 2 of 25 Veen et al. 2015Veen et al. , 2017, which is also sometimes referred to as control co-design (Garcia-Sanz 2019). Through integrated (topology) optimization, the dynamic behavior of the structure and the controller can both be adapted to accommodate each other in a more optimal manner, potentially resulting in a better closed-loop performance.
A large portion of existing research on integrated controller-structure optimization is focused on state-feedback controllers in the time domain, which determine their correction signals based on the state of the structure (e.g. positions, deformations, and/or velocities). The optimal controller in this case can be calculated algebraically as the minimizer of a linear quadratic control cost function, based on H 2 synthesis (a generalization of classic LQ/LQR/LQG theory) (Doyle et al. 1989;Anderson and Moore 1989). The result is an optimal controller balancing vibration levels and control effort over time. Most existing methods reformulate the integrated controller-structure optimization problem into a nested formulation, where an optimal controller is found algebraically during each structural design iteration (Haftka 1990;Fathy et al. 2001). The same linear quadratic cost function as the nested controller optimization can be used for the outer structural problem (Miller and Shim 1987). However, this approach is limited to truss problems with few design variables due to its significant computational effort, as the solutions are needed of an algebraic Ricatti equation and of additional Lyapunov equations for the gradients of each design variable. Alternatives in literature are based on minimizing combined strain and control energy in a steadystate setting (Ou and Kikuchi 1996a, b;Molter et al. 2013) or other (multi-)objective formulations (Zhu et al. 2002;da Silveira and Fonseca 2010). While these are computationally feasible for topology optimization, they do not directly relate to the integrated system performance. A more complete overview of different approaches is given in the review paper by Allison and Herber (2014).

Frequency domain control
In practice, state-feedback and linear quadratic optimal controllers in the time domain are rarely used for high-performance positioning systems. For positioning systems the tracking error, disturbance rejection, and noise attenuation are essential aspects to obtain high-precision. The quantification of these effects is difficult in the time domain (Doyle 1978;Zhou et al. 1996), but can be represented more easily in the frequency domain. This is one of the reasons that frequency-based proportional-integral-derivative (PID) controllers are the current industry standard (Munnig Schmidt et al. 2011).
In order to clearly describe the open challenges for integrated controller-structure optimization in the frequency domain, we will first discuss some aspects of classic control theory. The influence of disturbances and noise on the controlled structure is characterized by the sensitivity function S(j ) (Åström and Murray 2008), which is not to be confused with the design sensitivities. The sensitivity function is the transfer function between external disturbance d and output y (Fig. 1), which is dependent on the frequency . Here, the disturbances may, for instance, be external loads on the controlled system or motions of the measurement frame. The sensitivity function is defined as with the open-loop transfer function L(j ) = H(j )C(j ) of the controller and plant in series. The amplitude of the sensitivity function |S(j )| provides a bound on the disturbance rejection properties, of which a typical example is shown in Fig. 2a. Disturbances are attenuated by the feedback controller if |S(j )| < 0 dB , but they are amplified if |S(j )| > 0 dB . The controller is able to correct disturbances for frequencies below the bandwidth b , which ideally keeps the sensitivity function small at those frequencies and additionally ensures a small tracking error. However, the sensitivity function cannot be lowered indefinitely due to the waterbed effect: lowering the sensitivity function at certain frequencies leads to an increase at other frequencies (Munnig Schmidt et al. 2011). Therefore, peaks are to be avoided for frequencies above the bandwidth to prevent over-amplification of high-frequency noise. This is usually done by limiting the maximum value of |S(j )| to, for instance, 6 dB . Further details can be found in textbooks on control, e.g. Åström and Murray (2008); Munnig Schmidt et al. (2011).
Examples in the literature of integrated optimization for PID control are less common than optimization based on linear quadratic control. One example is the work by Albers and Ottnad (2010), who use a PID controller optimization nested within a structural topology optimization based on strain energy minimization. Here, load cases are iteratively updated based on the control action. However, this approach will not yield optimal performance, since the structure is optimized for a minimum strain energy instead of the integrated system performance.

Fig. 1
Controller and plant placed in a feedback loop, with the aim for the output y to track the reference signal r. The correction signal u generated by the controller is based on the measured error e. If tuned correctly, the closed-loop system is able to reject disturbances d A truly integrated approach is proposed by van der Veen et al. (2015van der Veen et al. ( , 2017, who optimize feedback controlled structures for a maximum bandwidth, subject to constraints on closed-loop stability and disturbance rejection. For the disturbance rejection, constraints are used that explicitly limit the sensitivity function |S(j )| below a certain threshold. Since the sensitivity function is a multi-modal function (Fig. 2a), a constraint is imposed on each individual peak value. The frequencies corresponding to the peak values cannot be calculated explicitly, so a numerical search algorithm must be used to locate the peak values (see, e.g. Bruinsma and Steinbuch 1990). Even though the peaks are found numerically and as long as the constraints are not dependent on the peak frequency, it is possible to calculate correct gradient information and use them as constraints in an optimization (Giesy and Lim 1993;Venini and Pingaro 2017;Delissen et al. 2020). A drawback of the approach of van der Veen et al. (2015van der Veen et al. ( , 2017 is that the number of peaks changes during the optimization, depending on the controller, the structure, and their interaction. Next to that, separate constraints need to be applied to ensure closed-loop system stability. As a result, integrated optimization including control requirements, such as closed-loop stability and disturbance rejection, remains an open challenge. For further insight into the behavior of the sensitivity function |S(j )| , an alternative interpretation is discussed. The sensitivity function can also be interpreted using the Nyquist curve of the open-loop transfer function L(j ) , as is shown in Fig. 2b. From Eq. 1 can be deduced that the reciprocal of the sensitivity function is equal to the distance from the Nyquist curve L(j ) to the critical point at −1 + 0j in the complex domain (from here on called the −1 point). A maximum in the sensitivity function therefore corresponds to a minimum distance between the open-loop transfer function L(j ) and the −1 point. This minimum distance is also commonly known as the modulus margin (Åström and Murray 2008).
Outside structural and integrated optimization, several techniques are available which focus on gradient-based tuning of controllers using the Nyquist curve, e.g. (Karimi and Galdos 2010;van Solingen et al. 2018). These apply geometric constraints on the Nyquist curve to enforce stability and disturbance rejection margins. An advantage of constraining the Nyquist curve is an enhanced flexibility in limiting both phase as well as amplitude of a transfer function, which is otherwise difficult to do. However, this approach is not suited for topology optimization, since the Nyquist curve L(j ) is sampled using a finite number of frequencies, where each sampled point has to be constrained in the complex domain. This easily results in thousands of constraints that each require a computationally costly (dynamic) finite element analysis.
Topology optimization has not been done yet using the Nyquist curve, although it may offer several advantages. There is a straightforward geometrical interpretation of disturbance rejection using the open-loop transfer function L(j ) as opposed to the sensitivity function |S(j )| , which is in closed-loop. Additionally, the closed-loop stability can be directly enforced by preventing encirclements of the −1 point (i.e. the Nyquist stability criterion in case of a stable open-loop system). This motivates the use of the Nyquist curve in controller design.

Robust formulation
An important requirement for practical design cases is the control on minimum feature size in the design. In topology optimization, this is generally done using a density filter in combination with a robust formulation (Bendsøe and Sigmund 2003;Wang et al. 2011). Erosion and dilation operations are performed on the design in order to generate multiple perturbed designs. By optimizing the design for worst case performance, it is made robust against uniform geometric deviations. The robust formulation indirectly ensures a minimum feature size in the design, dependent on the perturbation amount and the filter radius. Additionally, it helps in obtaining a binary void-solid design without intermediate densities, possibly also reducing the appearance of local eigenmodes (Pedersen 2000). A disadvantage is that the application of this method requires the solution of additional perturbed designs, which in the present setting each require the solution of a computationally costly eigenvalue problem. Furthermore, the integrated controllerstructure optimization as proposed by van der Veen et al. (2015van der Veen et al. ( , 2017 does not allow for aggregation of constraint values corresponding to the different perturbations. This is due to the fact that the number of peaks (and thus the number of constraints) may change due to the design perturbations, making it hard to determine which peaks to aggregate. The lack of aggregation results in the addition of many new constraints for each perturbed design, which all require the calculation of eigenmode design sensitivities. Thus, to apply the robust formulation to existing integrated optimization methods results in an unacceptable increase of computational effort by at least a factor three.

Approach and contributions
In this work, we present two main contributions towards integrated controller-structure optimization and application to more practical design cases: 1. Local approximation of the Nyquist curve using circles, which can be used in gradient-based optimization 2. An efficient robust method for dynamic topology optimization problems, that requires negligible additional computational effort These new methods are combined, tested, and demonstrated for the integrated controller-structure topology optimization of a motion system. Local approximation of the Nyquist curve In order to efficiently influence the shape of the Nyquist curve during optimization, local circular approximations are generated at each eigenfrequency. Using multiple circular approximations, the characteristic shape of the Nyquist curve is captured by simple geometric features. Finally, by geometric restriction of each circle in the complex domain, the global shape of the Nyquist curve can be influenced during optimization. This can be used to enforce closedloop stability (no encirclements around the −1 point) and robustness (minimum distance to the −1 point).
This method avoids the requirement of knowing the exact frequencies at peaks in the sensitivity function, or equivalently where the Nyquist curve of L(j ) is closest to the −1 point. Instead, locally approximated sections are used to describe the Nyquist curve close to the peak frequencies, which may be constrained away from the −1 point. Each (flexible) eigenmode in the mechanical model exhibits itself as a circle in the Nyquist curve (Fig. 2b), which together form its characteristic shape. This also prevents issues with a changing number of peaks in the sensitivity function, as the absence of a peak corresponds to a circle with a radius of zero at the corresponding eigenfrequency.
Circles in the Nyquist curve have historically been used to identify modal parameters of mechanical systems from experimental data (Kennedy and Pancu 1947;Miller 1978).
Here, the reverse process is exploited by fitting a circle related to each eigenfrequency in the dynamic system, using the corresponding modal parameters. Local approximation models are actively being researched in the field of control, where they are used to, for instance, approximate the H ∞ norm with limited experimental data (see, e.g, Tacx and Oomen (2021)). However, to our knowledge the current approximation-based approach proposed for integrated controller-structure optimization has not been studied before.
Computationally efficient robust formulation To apply the robust formulation to topology optimization with negligible additional computational effort, we propose to approximate both the eigenfrequencies and eigenmodes of the perturbed designs. This is critical, as these are both important to the closed-loop behavior of the system. Since the eroded and dilated designs are very similar to the nominal design, it may be assumed that their dynamic behavior is also very similar. After calculation of the eigenmodes of the nominal design, approximations of eigenfrequencies and eigenmodes for the perturbed designs are constructed from linear combinations of the nominal eigenmodes. This avoids solving additional eigenvalue problems for the perturbed designs.
The eigenmodes are used to construct a reduced-order model for each perturbed design. After this, closed-loop performance is evaluated for each of the reduced-order models, using the proposed local approximation method of the Nyquist curve. The fact that the number of circles does not change during iterations allows for aggregation of the constraints. Thus, the number of constraints is equal for an optimization with or without robust formulation, which prevents calculation of additional eigenmode design sensitivities.
In this work, the contributions are applied to integrated controller and topology optimization with focus on closedloop system stability and robustness margins on the disturbance rejection. The research is focused on mechanical linear time-invariant (LTI) and single-input single-output (SISO) systems, but many aspects can be generalized to a multi-input multi-output (MIMO) setting. A PID controller with a predefined structure is used for positioning of the mechanical system in a single motion direction. Although the controller-structure is known, its parameters are tuned during optimization. Fixed actuator and sensor locations are used in the system and modal damping is assumed. Furthermore, different problem variations are explored, such as optimization for position dependency using multiple sensor positions (i.e. single-input multiple-output (SIMO)) with the same controller and application of the proposed robust formulation.
The outline of this paper is as follows. First, in Sect. 2, the local approximation of the Nyquist curve using circles is explained and demonstrated on an analytical example. In Sect. 3, the topology optimization formulation is presented. Next, in Sect. 4, all modeling aspects are explained in detail, including the proposed efficient robust formulation for dynamic problems. In Sect. 5, the potential of the proposed methods is demonstrated using numerical examples. Finally, discussion and conclusions are given in Sects. 6 and 7, respectively.

Circle parametrization
In this section is explained how the circular-shaped local approximations for the flexible eigenmodes in the Nyquist curve L(j ) are constructed, based on theory from experimental modal analysis (Kennedy and Pancu 1947;Miller 1978). Given the transfer function in the Laplace domain L(s), the Nyquist curve is obtained with complex frequency s = j , which corresponds to a line along the imaginary axis. First, the general transfer function is given by its decomposition in first-order terms as with participation factors p i ∈ ℂ , and system poles i ∈ ℂ . It is possible to obtain this decomposition from any representation of the transfer function, for instance from a state-space model as is further explained in Sect. 4.3. For frequencies s = j in the proximity of a system pole i , the corresponding first-order term is assumed dominant, as its denominator becomes small. Therefore, where the local approximation L i (s) consists of a constant offset L i ∈ ℂ and a single first-order term. The offset L i contains the contributions of all remaining first-order terms at the frequency of interest and is calculated as ensuring interpolation of L(jIm i ) =L(jIm i ) . An illustration of a local approximation can be seen in Fig. 3.
From experimental modal analysis, it is known that a transfer function of the form in Eq. 3 results in a circle in the complex domain (Kennedy and Pancu 1947;Miller 1978). Its midpoint and radius are calculated, respectively, using An alternative and simplified proof of these relations, based on the theory of generalized circles (Schwerdtfeger 1979), is provided in Appendix 1. As the relations are all analytical, the derivatives can be calculated explicitly, for which the equations are given in Appendix 1. By constructing local circle approximations for each relevant system pole i , the important features of the Nyquist curve of L(s) can be described using simple geometry.
Note that the radius is non-differentiable when | | p i | | = 0 and additionally approaches infinity when Re i → 0 . The former case occurs when an eigenmode is not excited by the actuator or when it has no deformation at the sensor location. The latter case will normally not occur, since the real part of a complex pole is a finite negative value for a stable system Fig. 3 Illustration of one local circle approximation L (j ) (in blue) fitted to the Nyquist curve L(j ) in grey. The radius R, midpoint X, constant offset L , and the interpolatory point L(jIm( )) are indicated with damping. To account for non-differentiability when | | p i | | = 0 , and thus R i = 0 , a small perturbation is added as This ensures differentiability when R i = 0 by setting a (smooth) minimum radius of R min .

Analytical example
To demonstrate the principle of the local circle approximation, a double mass spring damper system is used, as is shown in Fig. 4a. To keep the equations simple, no controller is used for this example and circles are formed for the plant H(s) instead of the loop L(s). The transfer function H(s) for this system describes the relation between a force on either of the masses to a displacement on either. Mathematically, it is the superposition of two second-order systems (Gawronski 2004;Munnig Schmidt et al. 2011): Here, undamped eigenfrequencies are denoted Ω 1 , Ω 2 , relative damping factors 1 , 2 , and modal contributions 1 , 2 . These modal parameters can be calculated from the mass, stiffness and damping values of the double mass spring system in Fig. 4a and can be positive (e.g. collocated) or negative (e.g. non-collocated) (Gawronski 2004). This equation can be rewritten into a notation using system poles i and their conjugates i , becoming w h e r e t h e s yst e m p o l e s a r e c a l c u l a t e d a s case of an underdamped system. The Bode plot of this system in Fig. 4b shows the frequency-dependent amplitude and phase behavior of the transfer function H(j ) . Resonances can clearly be observed, which are located at the damped eigenfrequencies Im i in the Bode plot. Alternatively, the transfer function can be represented in the complex domain by a Nyquist plot, shown in Fig. 4c. Looking at the Nyquist plot, the circular shapes can clearly be identified, with their apexes with respect to the origin coinciding with the damped eigenfrequencies Im( ) . Note that the furthest point from the origin can be calculated explicitly, as opposed to the distance to the −1 point, which cannot be calculated.
Circles cannot be fit directly to the second-order systems, therefore the transfer function first needs to be decomposed into first-order terms. In this case, this is simply done by rewriting Eq. 8 into with corresponding participation factors From the four first-order systems in Eq. 9, only one firstorder system is assumed dominant around each of the excitation frequencies ≈ Im i and ≈ −Im i . Since only positive frequencies are of interest, the two approximated circles using Eqs. 3-5 can be described with radii and midpoints being with constant offsets Two circles are calculated and overlaid in Fig. 4c, demonstrating approximation with a close match to the original Nyquist curve. The approximation is most accurate for frequencies near the interpolation frequency Im i . From Eq. 11 can be seen that increasing the damping value or eigenfrequency would decrease circle radius, and the modal contributions i have a proportional effect on the radii. When eigenfrequencies are close to each other, they start interacting with each other. The mixing of modes results in non-circular shapes in the Nyquist curve. At this point the approximation in Eq. 3 is unable to fully capture the exact (10) behavior, as the influence of other modes ( H ) is no longer (close to) constant. Using the double second-order system, this is demonstrated in Fig. 5a for a system with positive modal contributions (e.g. in a collocated system). In this case the two second-order terms contribute in the same direction and the total response is larger than the approximations. The opposite happens when the signs of the modal contributions are opposed (e.g. in a non-collocated system), as can be seen in Fig. 5b, where the approximated circles are larger than the actual response. The two modes compensate each other and the total response becomes smaller. These examples are specifically tuned to show the effect of mode interaction, but these extreme cases may be rare to occur in an optimization setting. This will be investigated in Sect. 5 using numerical examples.

Constraining the Nyquist curve
Using the local circle approximations, the Nyquist curve can be parametrized in the complex domain using simple geometry. This is very useful for optimization problems where the Nyquist curve must geometrically be constrained in the complex domain. Distance to a point For instance, the closest distance h from a circle to a point ∈ ℂ is characterized as the distance to the midpoint of the circle minus its radius To calculate the distance furthest away, the radius is simply added instead of subtracted Distance to a line The shortest distance to a line characterized by unit normal direction n ∈ ℂ (with |n| = 1 ) and passing through the point can easily be calculated as Distance to an area By composing distances to lines and points, also the distance from a circle to an area can be characterized. For instance, the shortest distance to a wedgeshaped area bounded by two line sections with normals n 1 , n 2 ∈ ℂ intersecting in the point , is defined as where the conditions C 1 , C 2 , and C 3 indicate which of the three sections (lines or point) is closest to the position of X. These functions are mostly smooth and differentiable, except when point X coincides with or at an inflection point between two segments. However, these cases will numerically rarely occur, especially when these are used in constraints that serve to keep the point X away from .

Optimization formulation
In a closed-loop controlled system, the interplay between controller and structure determines the performance that can be achieved. The feedback system consists of a PID controller C(s) and the structure H(s), which contains a rigid body mode, placed in a loop, as is shown in Fig. 1.
From an optimization point of view, there is a trade-off between performance (bandwidth) and closed-loop stability. Stability can be determined by inspecting the closedloop poles, which must have a negative real part. Using the Nyquist stability criterion, closed-loop stability can also be interpreted with the Nyquist curve: for a stable closed As discussed in the Introduction, the modulus margin gives information on how close a system is to instability, and it also provides a bound on the influence of disturbances on the controlled structure. It is characterized as the closest distance of the Nyquist curve L(j ) to the −1 point.
To ensure a closed-loop system which is stable with a specified modulus margin, the locally approximated circles are used to constrain the trajectory of the Nyquist curve L(j ) . Constraints are defined to prevent the circle approximations, corresponding to the mechanical eigenmodes, from entering the wedge-shaped area offset by around the −1 point as indicated in Fig. 6. In this way, both stability and disturbance rejection are enforced simultaneously by the constraints.
Multiple Nyquist curves can be constrained in an optimization, for instance, to account for position-dependent dynamics. Many high-tech positioning systems consist of motion systems stacked in series to provide positioning freedom in additional movement directions or to achieve an extended range of motion. Since contactless measurements are often used (e.g, laser interferometry or eddy current sensors), a sensor fixed on a measurement frame therefore changes position relative to the measured object. As the measurement position affects the dynamics, it becomes position-dependent (van der Veen et al. 2017). In this work, the inclusion of multiple ( N out ) relative sensor positions is therefore also considered. This results in multiple SISO control loops, or rather SIMO, each requiring performance constraints on disturbance rejection.
For the topology optimization, a density-based formulation is used (Bendsøe and Sigmund 2003). The optimization formulation used in this work is stated as Fig. 6 The constrained area in the Nyquist diagram indicated in red, which ensures both closed-loop stability and a modulus margin of . The wedge-shaped area (Eq. 16) is shaded in dark red, which is offset with the modulus margin to the lighter shade of red where x represents the pseudo-density variables used for topology optimization and b the bandwidth, which is the tuning parameter of the PID controller. The number of flexible eigenmodes in the system is equal to N, for each of which circular approximations are constrained. The aim is to maximize the bandwidth b of the closed-loop system while keeping the volume V below a volume fraction v f of the maximum volume V max , and simultaneously ensuring the circle approximations remain outside of the wedge-shaped area using the distances h ij as defined in Eq. 16. These distances are related to the design variables through the radii R ij (x, b ) and midpoints X ij (x, b ) of the circular local approximations.

Optimization implementation and scaling
As optimizer, method of moving asymptotes (MMA) is used (Svanberg 1987). Constraint and objective scaling is critical to this method, so the original optimization formulation of Eq. 17 is reformulated as where the objective is normalized with the bandwidth at the initial iteration (0) b and a normalized design variable x is used to tune the controller. The constraints on the circles are scaled and normalized as Instead of directly using the bandwidth as a variable, it is exponentially scaled between the user-defined bounds [ min , max ] as This causes less sensitive behavior for parameter changes at low bandwidth, and makes all optimization variables ( x and x ) equally bounded between 0 and 1. The objective function is thus an explicit analytical function of the design variable x , hence the design sensitivity analysis is straightforward.
To ensure a feasible initial controller for a given initial structure x (0) (uniform densities equal to the volume fraction v f ), a separate controller optimization is performed prior to the integrated controller-structure optimization. The control variable x is found using the formulation of this preoptimization, given as which has its optimum at (0) b . This value is used as initial bandwidth for the integrated optimization given in Eq. 18.

Topology optimization parametrization
Since a density-based approach is used, the structural design variables x are first filtered using a standard density filter, resulting in the filtered design field x f (Bruns and Tortorelli 2001). The Young's modulus E i and density i of each finite element i in the domain E are obtained from the filtered design parameters using the following material interpolation The small minimum design density x min prevents the stiffness matrix from becoming exactly singular when design densities are zero.
The low mass-to-stiffness ratio in Eq. 22 for low densities largely prevents the occurrence of local eigenmodes (Olhoff and Du 2005). These are unwanted eigenmodes in low density areas, with low corresponding eigenfrequencies. Local modes are further prevented using a flood fill algorithm on the design vector x , removing any material that is disconnected or very loosely connected to actuator or sensor locations. Elements that are connected to the non-design domains through densities lower than 0.2 are recursively clipped to the maximum of their neighbors. In an extreme case, the disconnection of bodies results in additional rigid body modes at frequencies close to zero. These measures prevent undesired localized modes, improving the convergence of the optimization. , Next, the stiffness and mass matrices are assembled using the material properties E and . For this, a grid of bilinear quadrilateral finite elements is used, with a full integration scheme and a plane strain condition. The assembly is performed as where K 0 and M 0 represent element stiffness and (lumped) mass matrix, respectively, corresponding to unit material properties, and denotes the matrix assembly operation over the entire domain E.

Mechanical model
This section describes the different steps that are taken to calculate the circle radii and midpoints starting from the finite element matrices and the controller parameter. To guide the reader, a high-level visual guide of the different calculation steps is shown in Fig. 7. Furthermore, the equations to calculate gradients during each step of the optimization are given in Appendix 1.
From the mass and stiffness matrices M and K , a reducedorder model is first constructed using a truncated modal decomposition. This model approximates the displacement field u(t) by superposition of a number of eigenmodes i scaled over time with amplitudes q i (t) , denoted as where is the projection matrix containing all eigenmodes. The eigenfrequencies Ω i and corresponding eigenmodes i are obtained by solving the undamped eigenvalue problem for the lowest N + 1 modes using mass-normalization of the eigenmodes as T i M i = 1 . The lowest eigenfrequency Ω 0 = 0 rad/ms corresponds to the rigid body mode for the degree of freedom that is controlled by the PID controller.
The projection matrix is used to obtain the reduced equations of motion as where is a diagonal matrix containing the eigenfrequencies and the non-dimensional damping ratio. The unit input force vector, as exerted by the actuator, is denoted as f and the unit output displacement vectors as measured by the sensors with block-vector G , with N out columns for each sensor. The input is denoted u(t) and the outputs for all sensor positions y(t) , as indicated in Fig. 1. The transfer function of the plant becomes which describes the behavior between the input and N out outputs of the plant in the frequency domain.

Controller
A PID controller with additional low-pass filter is used for feedback control of the rigid body mode, which is defined by the control law Fig. 7 Overview of the modeling approach, starting from the controller design variable b and the finite element system matrices ( K and M ), which are dependent on the structural design variables ( x ). The flow diagram shows all steps taken in order to calculate the midpoint and radius of a circle approximation to the Nyquist curve with tuning parameters gain k and bandwidth b . This is a PID controller based on industry standard rules-of-thumb, with integral action until b ∕5 , phase lead between b ∕3 and 3 b , and first-order roll-off beyond 5 b (Munnig Schmidt et al. 2011;van der Veen et al. 2015). The Bode plot of this controller can be seen in Fig. 8. The bandwidth is used as a design variable during the optimization and the gain is calculated using Here, m is the mass of the system and the rigid body mode response of the plant is equal to

This ensures that the open-loop gain at the bandwidth
In current work, the gain value is chosen as k 0 = 1.1 , which ensures correct interaction between the controller and the rigid body mode (Munnig Schmidt et al. 2011). Note that the method is not limited to this specific control law and variations in control behavior and parametrization are possible.
The control law can be rewritten into state-space form where A c ∈ ℝ 3×3 , B c ∈ ℝ 3×1 , C c ∈ ℝ 1×3 represent the controller-structure in canonical form (also given in Appendix 1) (Skogestad and Postlethwaite 2001). The vector c contains the internal state of the PID controller and is of length 3. Note that other controller-structures might require a different number of internal states. The open-loop response is obtained by placing the controller (Eq. 31) and plant (Eq. 27) in series, connecting the output of the controller to the input of the plant. In the form of a state-space model in the time domain this becomes with system matrices (Fig. 7)

Modal decomposition
Before circles can be mapped to the open-loop transfer function L(j ) , the transfer function needs to be decomposed into first-order terms (Fig. 7). From the state-space model (Eq. 34), the poles can directly be obtained by an eigen-decomposition of the system matrix A as where matrix Q∈ ℂ M×M contains all eigenmodes of the (right) eigenvalue problem and matrix ∈ ℂ M×M has all the complex-valued poles i on its diagonal. Substitution into the transfer function of Eq. 34 yields (31) c(t) = A c c(t) + B c e(t), Here, the matrix P i denotes the participation factors of all input and output combinations for mode i. The participation factors can be calculated as for a general MIMO system, where the outputs are indexed with j and the inputs with k. The current application only considers one input, so the last index is omitted. The number of first-order terms equals the number of state variables M = 2N + 5 : two originating from each flexible eigenmode included in the reduced-order model, two from the rigid body mode, and three from the controller. The negative frequencies are not of interest and the poles corresponding to controller and the rigid body mode cannot be approximated by a circle. Therefore, N circles are fitted to the flexible modes and constrained in the complex domain. With this decomposition, the radii and midpoints of the circles can now be found using Eq. 5, which is the final outcome of the diagram in Fig. 7.

Efficient robust formulation
To apply the robust formulation in topology optimization, multiple perturbed designs are generated by erosion and dilation operations. This is done using the smooth Heaviside operator defined as with the filtered design field x f , resulting in a projected design x p (Wang et al. 2011). The parameter determines the edge contrast of the projection and the amount of dilation or erosion, where a value of = 0.5 corresponds to the nominal design. By choosing multiple different values of , multiple perturbed designs x p ( k ) are generated.
However, straightforward analysis of each of these perturbed designs (as described in preceding sections) results in an additional computational burden, as a reduced-order model has to be created for each design. This means that for each projected design the eigenvalue problem needs to be solved (Eq. 25), which is a very computationally intensive step in the analysis.
As an alternative, we propose to approximate the eigenfrequencies and eigenmodes of the perturbed designs, using the reduction basis with eigenmodes corresponding to the nominal design. This means that only the eigenvalue problem of the nominal design needs to be solved in each design iteration and it is assumed that the eigenmodes of the nominal model can be used to describe the behaviour of the other perturbed designs.
The eigenvalue approximation proceeds as follows: using the different perturbed designs x p ( k ) , instead of the filtered design x f , the corresponding mass M k and stiffness K k matrices are assembled using Eqs. 22 and 23. Next, the eigenvalue problem is solved (Eq. 25) using the mass and stiffness matrices corresponding to the nominal design, yielding the reduction basis . Instead of solving additional eigenvalue problems for the remaining projected designs, their system matrices are projected using the reduction basis belonging to the nominal design as Their dimensions correspond to the number of modes in the basis, thus K k ,M k ∈ ℝ (N+1)×(N+1) . These projected matrices are then diagonalized by solving the small eigenvalue problem resulting in The matrix k ∈ ℝ (N+1)×(N+1) is a diagonal matrix containing the approximate eigenfrequencies of the perturbed design, which are in fact Ritz values. The corresponding approximate eigenmodes are linear combinations of the nominal eigenmodes, calculated as V k . The system of equations for the perturbed designs now be found as The remainder of the analysis follows the same steps for each model as in Fig. 7, so first the controller is added to form the open-loop state-space model (Sect. 4.2). This is again decomposed into first-order systems by calculating the poles and participation factors (Sect. 4.3), after which circle approximations are formed for each eigenmode (Sect. 2.1). Finally, distances from the circles to the −1 point are calculated to form constraints. The effect of this approximation will be studied in Sect. 5.5.
After the calculation, each perturbed model k has different constraint values g ij,k for each of its circles, corresponding to mode i and sensor position j. To limit the number of constraints from the different models, they are aggregated using an induced aggregation function (Kennedy and Hicken 2015) for any constraint f = g ij . This function approximates the worst-case constraint value (i.e the maximum) between the perturbed projections, controlled by the parameter b. For a large parameter b this expression approaches the true maximum. This particular function is chosen because f = f k in case all values f k are equal. The robust parameter introduced in Eq. 38 is increased during the optimization, meaning that for the initial iterations all perturbed designs are similar, and so are their constraint values. With this choice of aggregation function, the constraint values are not underor over-estimated during the early phase of the optimization. The aggregation ensures the number of constraints does not increase when using the robust formulation, and thus no extra computational effort is required to calculate eigenmode design sensitivities (Appendix 1) (Lee 1999).

Case and settings
The numerical case that is used to demonstrate the method is shown in Fig. 9. To ensure a position-independent system, N out different sensor positions are defined at the measurement surface. Measuring at any of these locations and using that signal for feedback control should result in a closed-loop stable system with required disturbance rejection. In Table 1, the settings are listed as used in the optimization, where the material properties correspond to those of , aluminium. Furthermore, the maximum number of design iterations is limited to 200 to prevent excessive calculation times. For the optimization, MMA is used with default settings. A move limit of 0.05 is used on the design variables to prevent large steps and oscillations.

Sequential optimization
In order to compare performance, a reference case is presented based on a sequential optimization. First, the structure is found by maximization of eigenfrequencies, and subsequently the PID controller is optimized using the proposed method (Eq. 21). The optimization formulation used for the eigenfrequency maximization is given as in which g Ω is the objective function, defined as The superscripted variable g (0) Ω denotes the value at the initial design iteration. This formulation maximizes the harmonic mean of the first three eigenfrequencies, for which further details and sensitivity analysis can be found in, e.g. Ma et al. (1995); Delissen et al. (2022).
The resulting structure after optimization of the eigenfrequencies is shown in Fig. 10 and the subsequent controller optimization is able to achieve a bandwidth of 1.11 rad/ms . From the Nyquist plot in Fig. 10 can be seen . Fig. 9 The case used for the optimizations, which has one rigid body mode in vertical direction. The actuator forces (green arrows) act uniformly on the non-design domain (black) to the left. The top-right non-design domain represents the surface where accuracy is required.
Here, vertical displacements are measured at N out different sensor locations (blue arrows). The domain is of dimensions 300 × 100 mm , with an in-plane thickness of 300 mm , and is discretized into 210 × 70 elements that the controller satisfies closed-loop stability and disturbance rejection requirements. It can also be seen that only the sensor position at the tip (position 1) is limiting the bandwidth, of which the second eigenmode is touching the margin. Therefore, optimizing for different number of sensor positions will result in an equal bandwidth, provided the sensor location at the tip is included.

Integrated optimization
Using the proposed procedure, integrated optimizations are performed for different numbers of sensor positions N out = 1, 3, 6 . The designs and corresponding Nyquist plots are shown in Fig. 11. Mechanism-like structures can clearly be identified in the designs. The Nyquist plots show that all the designs meet the requirements on closed-loop stability and modulus margin. However, not all designs contain binary zero-and-one densities that can directly be interpreted. Especially the design for one sensor position ( N out = 1 ) contains large areas with intermediate densities.
The design with six sensor positions also contains some areas with intermediate densities. This might be functionally interpreted as a 'rubber band' with a specific stiffness to tune the system dynamics. An overview of the achieved performance, as compared to the sequentially optimized design, is shown in Table 2. All the designs optimized with the integrated approach have a bandwidth about a factor 3.5 higher than the design optimized for eigenfrequencies. Moreover, the eigenfrequencies are significantly lower for the integrated optimizations, which clearly demonstrates that the system with the highest bandwidth does not necessarily need maximized eigenfrequencies.
The integrated approach is able to achieve a high bandwidth, relatively close to the eigenfrequencies. This can be explained using the Bode diagrams in Fig. 12, in which the dynamic response of the design optimized for integrated performance with N out = 3 is shown. Fig. 12a shows the dynamic response of the plant. The main difference between the sequential and integrated optimized design is the fact that the integrated design generally has smaller resonance amplitudes. This allows a controller with higher bandwidth, since the eigenmode is less dominant. The effect of a controller with a higher bandwidth can be seen in the Bode plot of the open-loop system (Fig. 12b). A controller with a higher bandwidth adds more gain to the system, hence the amplitude of the integrated design is higher than that of the sequential design, particularly in the low frequency range below the bandwidth. At higher frequencies, the amplitudes of the peak frequencies are about the same height in the open-loop responses, which is a result of the disturbance rejection constraints. The first eigenmode of the integrated design creates a resonance peak with a very small amplitude (Fig. 12a), around 5.2 rad/ms . This small amplitude means that the actuator is unable to 'affect' this mode and/or it cannot be 'seen' by the sensor (i.e. uncontrollable and/or unobservable). For the controller it seems as if this mode does not exist, therefore it is not limiting bandwidth. Inspecting the mode shapes of the integrated designs in Fig. 13, this effect can clearly be seen. For some modes, the actuator is virtually at a standstill, meaning that the mode is not excited by the actuator. For other modes, the location corresponding to the sensor is at a standstill, which means the sensor does not measure the mode.
The convergence history of the three designs is shown in Fig. 14. Especially the design for N out = 1 shows significant oscillations. In the Nyquist curves of subsequent design iterations shown in Fig. 15, the circle corresponding to the first eigenfrequency flips its direction. This flipping is caused by the actuator or sensor displacement crossing zero and changing sign. Since the modes have a very small excitation amplitude (Fig. 13), the controller is able to attain a very high gain. A small variation in the design then causes a small change in the mode shape, which eventually has a large effect on the system, due to the high control gain. The designs for 3 and 6 sensor positions exhibit less oscillations and a smoother convergence. Due to the addition of multiple sensor locations, the complexity of the optimization problem is increased and involves more trade-offs, leading to designs which are less sensitive to small variations.

Comparison with explicit peak constraint
To demonstrate the added value of proposed method, also an optimization based on the method of van der Veen et al. (2015,2017) is implemented. First, the frequencies corresponding to peaks in the sensitivity function are located numerically in each design iteration, after which they are used as constraints. Additionally, an explicit constraint ensuring closed-loop stability must be added, to prevent the Nyquist curve from encircling the −1 point. This is done by limiting the (smooth) maximum of all real parts of the closed-loop poles below zero, thus ensuring all poles are in the left half plane. We refer to the original publications for the full description, as our interest here is primarily in the comparison with our proposed approach.  Fig. 16. Although a structure may be recognized in the designs, the structural features look rather irregular and contain substantial amounts of intermediate densities. The design considering one actuator position faces similar convergence issues as the proposed integrated optimization (Fig. 11a), where small changes in design and mode shape cause oscillations. Next to that, the optimization for multiple sensor positions results in infeasible designs. During the design iterations, one of the Nyquist curves loops around the −1 point, indicating closed-loop instability (Fig. 16b, c). In this situation, two constraints are conflicting: the stability constraint requires the Nyquist curve to pass on the right side of the −1 point, but the peak constraint prevents this by requiring the curve to stay outside the circular margin. A change of design or bandwidth will thus violate at least one of the constraints, making it difficult to escape this situation. The proposed method does not face these issues, as stability is ensured implicitly

Robust formulation
For the application of the robust formulation to the proposed integrated optimization, the effect of the robust parameter is studied first. This parameter controls the amount of dilation or erosion of the design. Both the eigenfrequencies and the distances from the circle approximations to the −1 point change as a function of , as is shown in Fig. 17. Here, a distinction is made between the approximate responses using the nominal eigenmodes, as explained in Sect. 4.4, and those evaluated exactly, using eigenmodes corresponding to the perturbed designs. As can be expected, the error between the exact and approximated responses deviates more as the design is perturbed further away from the nominal design at = 0.5. Another observation that can be made in Fig. 17, is the fact that the distances to the −1 point are not monotonically increasing or decreasing, as is the case for compliance problems (Sigmund 2009). The lack of a monotonic behavior means that the worst-case design is not necessarily coinciding with extreme values of . In this work, a value of = 0.5 ± 0.05 is used, for which it can be assumed that the worst-case performance is likely to be included by evaluating three designs, at = 0.45 , 0.5, and 0.55. However, for larger perturbations and given the non-monotonic behavior, evaluation of the worst-case performance might require more than three designs to be analyzed.
For the robust optimization procedure, the edge contrast parameter is gradually increased from 1.0 to 20.0 during design iterations 50-180. Also, a filter radius of 8 elements is used to ensure a large minimum feature size. The resulting designs of the robust integrated optimizations, using the proposed method, are shown in Fig. 18. All designs have clear boundaries between void and material, which is a characteristic property of the robust formulation. At some locations, hinges appear with intermediate densities to provide a low-stiffness connection. As for their performance, the design with one sensor position achieves a bandwidth of b = 3.5 rad/ms , and the two other designs a bandwidth of 1.9 rad/ms . A design robustly optimized for maximum eigenfrequencies, with similar settings, is found to achieve a bandwidth of 1.1 rad/ms . This means the performance increase of robust integrated optimization is still significant.
Although the bandwidth of the design with one sensor position is very high, its disturbance rejection requirement is not met, as can be seen in Fig. 19a. This can again be explained by the dynamic response being very sensitive to (small) design variations. For the designs optimized with 3 ( Fig. 19b) and 6 sensor locations, the disturbance rejection requirements are satisfied for all three design perturbations and at all sensor locations.
The convergence properties are also improved using the robust formulation, as is seen in Fig. 20. Small oscillations are still present, but significantly less than without the robust formulation (Fig. 14). A lower final bandwidth is attained for all designs, compared to the results from optimizations without the robust formulation. However, this is counterbalanced with an increased robustness against geometric perturbations and the added control on minimum feature size.
The time required for this optimization using the approximated robust formulation is 55 min. The same optimization without robust formulation requires 49 min in total. This is for a total of 200 design iterations and 3 sensor positions, on a standard corporate laptop with an Intel Core i7-6600 processor. The difference of 6 min (12% extra calculation time) is required for the Heaviside projections, approximation of perturbed eigenmodes, and analysis of the reducedorder model for the local approximations. This is significantly less than the factor 3 that would be required without the proposed method approximating the eigenmodes for the perturbed designs.

Discussion
In the numerical examples, no issues with severe mode interaction and ill-fitting local approximations are encountered during the optimizations, which were initially identified as potentially harmful for the optimization (Fig. 5). The current work is focused on monolithic structures made out of metal, with a low damping coefficient . Mode interaction might increase for applications with higher damping, which might lead to issues in convergence. More research is required into the effect of damping and mode interaction on the optimization.
In some of the designs from the numerical examples (Fig. 11), intermediate density values are present. This can make the interpretation into solid and void material difficult from a practical engineering point of view. Partly, the intermediate densities might be caused by the fact that material penalization is harder for the integrated optimization problem than, for instance, compliance minimization or eigenfrequency maximization problems. For the latter, penalization is straightforward to apply, as more mass means a lower eigenfrequency and less stiffness means a higher eigenfrequency. For the integrated optimization problem it is difficult to make this distinction, as frequencies and modes are specifically tuned with respect to the controller behavior. As can be expected, application of the robust formulation helps in reducing the amount of intermediate densities in the designs (Fig. 18). Nonetheless, the reduction of intermediate densities for complex topology optimization problems remains an important topic for research.
In the examples, a discrete number of sensor positions is used to account for position-dependent dynamics. However, position dependency occurs over a continuous line or even a surface in reality. The current implementation ensures accuracy at several locations, but it does not account for any other locations also requiring accurate positioning. Even when locations may not be sensed, their accuracy may still be important to consider in the design of a motion system. Moreover, the inclusion of an increasing number of sensor positions leads to a growth in the number of constraints for the optimization. As eigenmode design sensitivities are required for each constraint, this directly increases computational cost. Therefore, future research should focus on methods to include many sensor locations without excessive computational cost.
The presented framework potentially allows for frequency-dependent constraints. This could be used, for instance, to incorporate more stringent limitations on highfrequency eigenmodes. Another possibility is to include the influence of time delay in the control loop, which has a frequency-dependent effect on the phase. Following the same approach based on the Nyquist curve could allow for consideration of this effect.
In reality, many systems are MIMO, for which the methods presented in the current work may also be used (e.g. sequential loop closing (Skogestad and Postlethwaite 2001)). However, for a coupled MIMO system multiple loops interact with each other, which requires further analysis of closed-loop behavior. This is done in the method of van der Veen et al. (2017), but it requires additional constraints and thus significantly more computational time. Efficient extension to MIMO systems therefore remains an open issue for future research.

Conclusion
A novel approach to integrated controller-structure topology optimization is proposed, with the aim of optimizing the structural design and tuning the controller parameters simultaneously for closed-loop system performance. We introduce a flexible framework enabling local approximation of the Nyquist curve using circles. These allow simple formulation of geometric constraints on the Nyquist curve in the complex domain, making them suitable for gradient-based optimization. In this manuscript, the approximating circles are used to constrain the shape of the Nyquist curve in an integrated controller-structure optimization, which enforces stability and disturbance rejection properties of the closedloop system. The approach is general and can be extended to other control objectives or constraints that can be expressed by the Nyquist curve.
From the numerical examples, it can be seen that the proposed method is able to greatly improve closed-loop system performance. For the studied problem, the state-of-the-art method in literature (van der Veen et al. 2015) is not able to converge to feasible designs due to conflicting stability and disturbance rejection constraints. In the proposed method, this problem does not occur as stability and disturbance rejection are ensured simultaneously by geometrical restriction of the local approximations of the Nyquist curve. Using numerical examples, the integrated optimization achieves improvements up to 350% in terms of bandwidth compared to sequential optimization, while ensuring a sufficient modulus margin.
Also position-dependent dynamics is considered, by the addition of constraints on multiple SISO Nyquist curves corresponding to different sensor positions. Not only does this lead to a structure and controller that can be used at each of the sensor locations, it also improves convergence properties of the optimization. Optimizing for only one sensor position results in designs that are very sensitive to small variations in mode shape around the actuator and sensor positions, leading to severe oscillations. By optimizing for multiple sensor locations, this detrimental effect is noticeably reduced.
Furthermore, a computationally efficient robust formulation is introduced, approximating the dynamics of the eroded and dilated designs. It allows for analysis and optimization of perturbed designs without significant additional computational cost, instead of a threefold increase using the conventional approach. The validity of the approximation is demonstrated for a range of design perturbations. For larger design perturbations extra care is required, because non-monotonic behavior is observed for the modulus margins, potentially resulting in interior worst-cases for robust optimization. Using the proposed formulation, the obtained designs are more robust against geometric deviations, a length scale is imposed, and a positive effect on optimization convergence is observed. Furthermore, the approximation-based robust formulation is not limited to the current application, but may also be used for other types of topology optimization involving dynamics.
There are several gaps to bridge in order to arrive at realworld systems with ultimate performance, such as incorporation of time delay, MIMO control, ensuring accuracy over large surfaces, and further reduction of computational effort. Despite the remaining challenges, this work provides a step forward in computational design methods for next generation high-precision motion systems.
Next, the variable L (Eq. 4) is directly substituted into Eq. 5. This leaves the sensitivities of the interpolation point L = L(jIm i ) ∈ ℂ , which are simply Finally, the partial sensitivities of the system pole i ∈ ℂ are Note that L also depends on i , for which a second contribution is calculated below in Eq. 70.

Open-loop response
From here on, we will switch to the matrix form of the sensitivities to calculate the contributions of different inputs and outputs of the system simultaneously. The matrix L(jIm i ) = L ∈ ℂ N out ×N in contains the various open-loop responses from different inputs to outputs.
Using the solution U and adjoint solution Y , calculated as the sensitivities can be calculated. First of the system matrix A as in which the sensitivity with respect to a matrix is a matrix whose entries are h . The modal decomposition basis Q also depends on the system matrix A , so the second contribution to this sensitivity is given below in Eq. 77. Secondly, the sensitivities of the input matrix B and output matrix C are (64) Besides L , the participation factors P i also depend on B and C , for which the sensitivity contributions are given below in Eq. 72. Finally, the sensitivity of the system pole i (the second contribution, next to the one given in Eq. 65) is calculated as in which the operation A ∶ B = ∑ i ∑ j A ij B ij is a double contraction.

Appendix B.4: Participation factor
Participation factor (Q, B, C) → P i , Eq. 37 To simplify the notation of sensitivities, the following notation is used here: Using this notation, Eq. 37 can be rewritten as P i = Cq iq T i B. The sensitivity contributions (the second contribution besides Eq. 68) of the input matrix B and output matrix C can be calculated using Here, the outer product is defined as A = u ⊗ v ↔ A ij = u i v j .
The sensitivities with respect to mode q i and q i (which is in fact the left eigenmode) are calculated as Since q i and q i are related to each other through Q −1 , the sensitivities with respect to the matrix Q can now be calculated as This is a combination of the adjoint for ĥ q i and then adding h q i in the correct column (i).

Modal decomposition
A → ( i , Q) , Eq. 35 Following the approach of Lee (1999), the eigenvector and eigenvalue sensitivities can be calculated using the adjoint method. Here again the matrix contains the sensitivities of the eigenmodes in its columns.
First the adjoint saddle-point problem needs to be solved for each modal pair ( i , q i ) , which is denoted This needs to be done for all eigenmodes with h q i ≠ 0 and/ or h i ≠ 0 . Now the sensitivity with respect to the system matrix is calculated as This is the second contribution, besides the one given in Eq. 67. Note that these sensitivity calculations are only valid with eigenmode normalization as q T q = 1.

Open-loop state space
( , , A c , B c , C c ) → (A, B, C) , Eq. 33 The matrices h A , h B , and h C are divided into partitions such that they match the division in Eq. 33: The row and column sizes in this division are equal to (3, N + 1, N + 1) . Using this division, the sensitivities with respect to the input variables can be calculated as

Modal analysis
(K, M) → ( , ) , Eq. 25 The modal sensitivities are calculated again using the adjoint method of Lee (1999), but now for a generalized eigenvalue problem with symmetric system matrices ( K and M ). First the adjoint saddle-point system needs to be solved for each modal pair (Ω i , i ) , which is a time consuming step as the size of this problem is (N dof + 1) × (N dof + 1) , where N dof is the number of unknowns in the finite element analysis. The adjoint problem is given as To represent the control law in Eq. 28 as a state-space model, the controllable canonical form is used, which results in the following matrices Now the sensitivities with respect to the control parameters can simply be calculated as Further sensitivity calculations (e.g. the dependency of k on b , Eq. 29) are straightforward.

Appendix B.8: Robust formulation
For the robust formulation, we start with the sensitivities h W k (with W k = V k ) and h k , which are equal to, respectively, h (Eq. 82) and h (Eq. 83) that were derived earlier.
The sensitivities of the linear combination factors ( V k ) used to combine the different eigenmodes are calculated as Yet another eigenvector adjoint system is constructed for the robust modal decomposition of Eq. 41 as Here, v i denotes the ith eigenvector of the system, which are normalized as v T iM v i = 1 . The different eigenvectors are the columns of V k , thus Using the solution pairs ( i , i ) the sensitivities with respect to the reduced-order matrices K and M can be calculated as Both the projection operation (Eq. 39) and the linear combination of eigenvectors ( W k = V k ) contribute to the sensitivities with respect to with Now, finally the large eigenvector system is solved again (Eq. 84) and the sensitivities with respect to the design variables x j are calculated as