Semi-active $\mathcal{H}_{\infty}$ damping optimization by adaptive interpolation

In this work we consider the problem of semi-active damping optimization of mechanical systems with fixed damper positions. Our goal is to compute a damping that is locally optimal with respect to the $\mathcal{H}_\infty$-norm of the transfer function from the exogenous inputs to the performance outputs. We make use of a new greedy method for computing the $\mathcal{H}_\infty$-norm of a transfer function based on rational interpolation. In this paper, this approach is adapted to parameter-dependent transfer functions. The interpolation leads to parametric reduced-order models that can be optimized more efficiently. At the optimizers we then take new interpolation points to refine the reduced-order model and to obtain updated optimizers. In our numerical examples we show that this approach normally converges fast and thus can highly accelerate the optimization procedure. Another contribution of this work are heuristics for choosing initial interpolation points.


Introduction
Consider a vibrational system described by a system of a second-order differential equations where M, K ∈ R n×n are the symmetric and positive definite mass and stiffness matrices, respectively.The internal damping C int ∈ R n×n is usually taken to be a small multiple of the critical damping denoted by C crit , that is see, e. g., [7,8,32].The presence of the positive definite internal damping C int ensures that the homogeneous system (1a) is asymptotically stable.Moreover, we have B 2 ∈ R n×p , E 2 ∈ R n×m , C 2 ∈ R p×n , and H 1 ∈ R ×n .The vector q(t) ∈ R n represents the displacements of the masses, u(t) ∈ R p is the control input, and y(t) ∈ R p is the measured output.Moreover, w(t) ∈ R m is an exogenous input signal, while z(t) ∈ R denotes the performance output.In principal, z(t) can also include an additional part which corresponds to the velocities of the masses q(t), but in this paper only the states are of interest.
In this paper we assume negative linear feedback corresponding to a linear damper of the form where G = diag(g 1 , g 2 , . . ., g p ) ∈ R p×p is the diagonal damping matrix with non-negative parameters encoded in a vector g = g 1 g 2 . . .g p T .The entries g i , i = 1, 2, . . ., p represent the friction coefficients of the corresponding dampers, usually called gains.In this paper we will consider the gains to be scalar variables which need to be optimized.
Here we assume that we have collocated inputs and outputs, i. e., C 2 = B T 2 .By using the feedback control law as in (3) we obtain the closed-loop system (For simplicity, we omit the dependency of q(•) on g and only add it to z(•).) M q(t) + C(g) q(t) + Kq(t) = E 2 w(t), (4a) where C(g) := C int + B 2 diag(g 1 , g 2 , . . ., g p )B T 2 is symmetric and positive definite.This implies that the unforced closed-loop system is also asymptotically stable.More details regarding system stability and model description can be found in [5,10,8].
With the substitutions x 1 (t) := q(t), x 2 (t) := q(t) and x(t) := x1(t) x2(t) we obtain a first-order representation of the closed-loop system Using the Laplace transform we obtain the closed-loop transfer function of (4), which is given by which is a real-rational matrix-valued function for each g ∈ (R ≥0 ) p , where R ≥0 denotes the set of nonnegative real numbers.The damping optimization problem has been investigated widely in the literature, for example in the books [17,33,26] where one can also find overview of different damping criteria.But even nowadays there are lot of open problems, especially from the computational point of view.One case studied in the past is the passive damping.This means that for given positive definite matrices M and K, the "best" damping matrix C(g) (with a certain structure) should be determined, such that the solution trajectories of the stationary system M q(t) + C(g) q(t) + Kq(t) = 0 have an "optimal" transient behavior.In this setting, gain optimization can be a computationally very demanding task, especially if the damper parameter values should be optimized for many different damper locations.There exist a number of different methods such as [16,12].Moreover, there are methods based on dimension reduction such as [7,8,30,31], where the minimization of the total average energy is the objective.
In this work we consider the semi-active damping optimization problem: we want to determine "the best" damping matrix C(g) which will minimize the output z(g, •) under the influence of the input w(•).The influence of the input on the output can be measured in several ways.There exist several optimization criteria such as the impulse response energy, the peak-to-peak performance, and the energy-to-energy performance.An overview of the different criteria can be found in [10].
Damping optimization using the impulse response energy leads to the minimization of the H 2 -norm of the closed-loop transfer function.This requires to solve an associated Lyapunov equation numerous times, which can make the optimization process computationally inefficient.In order to accelerate it, in [5] the parametric (subspace accelerated) dominant pole algorithm is used for the approximation of the impulse response energy.Moreover, in [29] the authors propose an efficient optimization approach using structure-preserving parametric model reduction based on the iterative rational Krylov algorithm (IRKA).There, several adaptive sampling strategies are used to obtain good approximations with respect to the H 2 -norm, aligning well with the underlying design objectives.
The method used in this paper is also based on adaptive sampling, but it differs from the approach for the H 2 -norm case.However, here we consider the energy-to-energy performance of the closed-loop system which is defined by where (w, x, z(g, •)) are measurable and square-integrable solution trajectories of the closed-loop system with parameters g.Hence, the energy-to-energy performance can be interpreted as the worst-case amplification of the energy of an exogenous input signal in the performance output.It is well-known (see, e. g., [34,Chap. 3]) that this criterion is equivalent to the H ∞ -norm of the closed-loop transfer function, that is Since the closed-loop system is asymptotically stable, then for all g ∈ (R ≥0 ) p , the functions F (g, •) are in the space where This space is equipped with the H ∞ -norm see [34,Chap. 3].Thus, in the setting of H ∞ damping optimization the problem is to determine optimal gains (under the assumption that the minimum exists).The optimization of damping based on the closed-loop H ∞ -norm can be a computationally challenging task.First, the optimization problem we consider here is nonlinear, nonconvex, and nonsmooth.The latter means that g * may be attained at a point, where the objective function g → F (g, •) H∞ is not differentiable.This problem occurs frequently in fixed-order H ∞ -control, see [19].Specialized methods for such optimization problems are available [22,23,15].These are modifications of quasi-Newton methods but in case of a nonsmooth optimizer, the H ∞ -norm (and its gradient) may have to be evaluated for a lot of parameter values.
Moreover, the calculation of the H ∞ -norm for a fixed parameter can be expensive, especially if the state-space dimension of the closed-loop system is large.The latter problem has been addressed by various works [18,9,25].However, in this paper we will modify the idea from [2] to our problem.In the latter work, the transfer function is approximated by interpolatory reduced-order models.Then the projection spaces are updated using the information of the point on the imaginary axis, at which the maximum of the largest singular value for the reduced order model is attained.This leads to an iterative algorithm that can be shown to have a superlinear rate of convergence to a local maximizer.In [3], this approach has been extended to the problem of minimizing the H ∞ -norm of a parameterdependent transfer function which is particularly useful in our context.With this approach, in each step we construct a parametric reduced-order model with reduced transfer function F (g, s) and compute Hence, in each step, one optimization of the L ∞ -norm (for the definition see ( 9)) of a parametric reduced transfer function and one large-scale non-parametric H ∞ -norm computation is needed.Then we choose g, i ω as a new interpolation point to update the reduced-order model and repeat this process until convergence.
However, in the context of damping optimization, the choice of the interpolation points is particularly subtle.Therefore, the main contribution of this work consists of deriving two heuristics for choosing initial interpolation points.We further do not only demonstrate the effectiveness, but also the numerical challenges of these ideas on various numerical examples.

Preliminaries
In this section we discuss some preliminary results for our method.In particular, we give some more details on how to compute the H ∞ -norm since this is the core idea of our algorithm.This method is derived in detail in [2] and we closely follow the presentation in [1].Moreover, we derive an analytic formula for the gradient of the H ∞ -norm with respect to the parameters g.

Computation of the H
is given.To compute F H∞ , in [2] reduced functions with and projection matrices V, W ∈ C n×k where k n are iteratively constructed.Since s E − A is of very small dimension compared to sE − A, the L ∞ -norm of F can be easily obtained.Assume that F is a transfer function that is bounded on the imaginary axis.Then the L ∞ -norm of F is defined by Note that, if F ∈ H m× ∞ , then we have F H∞ = F L∞ .In the method a sequence of such reduced functions F 1 , F 2 , . . . is constructed, whose L ∞ -norms converge to the H ∞ -norm of the original transfer function F .To obtain the reduced functions, the interpolation technique of [4] is employed.If m = and r interpolation points iω 1 , . . ., iω r are given, then F r is defined by F as in (7) and (8) (10) These choices for V r and W r imply that the Hermite interpolation properties are satisfied.It is important that V r and W r have the same number of columns so that the inverse s E − A −1 exists.This condition is usually violated by the projection matrices in (10) if m = .By using the projection matrices instead, this problem is solved and the Hermite interpolation properties (11) are preserved, see [2] for the derivation.Alternatively, a regularization procedure may be carried out to achieve this goal.Now, since the transfer function F r is constructed by matrices of small dimension, its L ∞ -norm can be obtained by well-established methods for the small and dense case which are reported in [13,11,6].With the point iω r+1 at which the L ∞ -norm of the current iterate F r is attained, the projection matrices are updated as In case of convergence, using the Hermite interpolation conditions (11), a superlinear convergence rate to a local maximizer of σ max (F (i•)) can be shown, see [2] for the details.A MATLAB implementation of this procedure, called linorm subsp, is publicly available * .

Gradient of the H ∞ -Norm
For the derivation of the gradient of the H ∞ -norm we make use of the following lemmas: Let σ(t) be a singular value of L(t) converging to a simple nonzero singular value σ 0 of L 0 as t → 0.Then, σ(t) is analytic near t = 0 and dσ(t) where u 0 and v 0 with u 0 2 = v 0 2 = 1 are, respectively, the right and left singular vectors of L 0 corresponding to σ 0 .
[24] Let s 0 ∈ C not be a pole of the transfer function can be expanded into a Laurent series at s 0 as From the two lemmas above we obtain the following result.
[9] Let s 0 = γ 0 + iδ 0 ∈ C not be a pole of the transfer function F (•). Furthermore, assume that the largest singular value of F (s 0 ) is simple with associated right and left singular vectors u 0 and v 0 satisfying u 0 2 = u 0 2 = 1 and define h(γ, δ) Then the gradient of h with respect to the variable γ is given by Now we derive a gradient formula for the function g → F (g, •) H∞ .For the derivation we make the following assumptions.Assumption 4. Assume that for some fixed g 0 , we have A2.The maximum singular value of F (g 0 , iω 0 ) is simple with associated right and left singular vectors u 0 and v 0 satisfying u 0 2 = v 0 2 = 1.
Assumption A1 ensures that the optimal frequency ω 0 at which the H ∞ -norm is attained is also uniquely determined in a neighborhood of g 0 .Therefore, in a neighborhood of g 0 there are no jumps in the points at which the H ∞ -norm is attained and thus there are no "kinks" in the graph of the function g → F (g, •) H∞ .Assumption A2 ensures that we can apply Lemma 3 to compute the gradient of the maximum singular value function using the corresponding singular vectors.Both conditions together imply that the function g → F (g, •) H∞ is differentiable in a neighborhood of g 0 and that we can make of use of formula (12) to compute the gradient of the H ∞ -norm as outlined below.
For g = g 1 g 2 . . .g p T and all j = 1, . . ., p, we define the matrices G j = diag(g 1 , . . ., g j−1 , 0, g j+1 , . . ., g p ) and Then, we can rewrite F (g, s) as for some ω 0 ∈ R ≥0 and let the corresponding right and left singular vectors be u 0 ∈ C and v 0 ∈ C m .Let Assumption 4 hold accordingly for g 0 = g.Then with Lemma 3 we obtain 3 Our Approach

Algorithm Description
In this subsection we present our adaptive interpolation approach for damping optimization which has been adapted from [3].The basic algorithm is given in Algorithm 1.Note that we have formulated the algorithm only for the case m = to simplify its formulation.If m = , then the updates of the projection matrices must be adapted as in Subsection 2.1.Some words of explanation are in order.To initialize the algorithm, in lines 1-4 initial projection matrices are constructed.This is done by sampling the parameter domain (R ≥0 ) p and the imaginary axis appropriately.Since the choice of good sampling points can be a difficult issue, we will give more details about this in Subsection 3.2.
In lines 5-11, the optimization loop is performed.In line 6, we construct a reduced-order model and its associated transfer function using the projection matrices.Since the reduced-order model is of low dimension, its parameters can be efficiently optimized in line 7.This is a nonconvex and nonsmooth optimization problem.The nonsmoothness arises from the fact that it might well happen, that the function g → F (g, •) H∞ is continuous but not differentiable for g = g * .This is the case, when for F (g * , •) there exist multiple frequencies at which the H ∞ -norm is attained (i.e., Assumption A1 is violated).Luckily, there exist optimization algorithms that are specialized for this kind of problem.Therefore, we make use of the method presented in [15] which can handle such problems and which has been implemented in the MATLAB package GRANSO † .This software only requires a function handle for the evaluation of the mapping g → F (g, •) H∞ and its gradient as discussed in Section 2.Moreover, constraints such as the elementwise nonnegativity of g can be easily incorporated into GRANSO.Then in line 8 we evaluate the H ∞ -norm of the original transfer function at the computed optimal parameter value to find a new optimal frequency for sampling and updating the reduced-order model.This is done using the subspace method of [2] for computing the H ∞ -norm of a large-scale system which is implemented in the MATLAB code linorm subsp ‡ .Finally, in line 9 of Algorithm 1 we update the projection matrices according to the optimal parameters we have just obtained.Note that our projection Algorithm 1 Greedy algorithm for semi-active H ∞ damping optimization Input: The matrices M, K, C int , B 2 , E 2 , H 1 defining the system (1) and the associated matrices B, C and the function D(•, •) as in (5), some initial parameters g (1) , g (2) , . . ., g (ν) ∈ R p Output: Optimal gains g * = arg min g∈(R ≥0 ) p F (g, •) H∞ , where F is the closed-loop transfer function in (5).

9:
Update the subspaces and set Orthonormalize the columns of V j and W j , respectively.11: end for 12: Set g * := g (r) .spaces are constructed in a way such that for i = 1, . . ., ν, j = 1, . . ., r, and = 1, . . ., p we fulfill the Hermite interpolation conditions This also directly leads to Hermite interpolation conditions between the functions σ max (F (•, i•)) and σ max ( F r •, i•) .In particular, we have This property is useful to show a superlinear rate of convergence [3], at least in the case of convergence to a differentiable minimizer.Finally, for avoiding ill-conditioning of the reduced problems, we orthonomalize the columns of V j and W j , respectively in line 10.
In this work we consider systems with a larger number of inputs and outputs.In order to avoid a fast growth of the subspace dimensions in our projection approach, we make use of tangential interpolation.That is, with v 1 ∈ C , w 1 ∈ C m being the right and left singular vectors associated with the largest singular value of F ( g (j) , i ω j ) we replace Step 9 in Algorithm 1 by with the function D(•, •) as in (5).This is similarly done in the projection approach for the computing the H ∞ -norm.In this way, the subspace dimension grows by one in each step which makes these computations feasible even for larger iteration counts.On the other hand, the convergence analyses of [2,3] are not valid any longer, since we may loose the Hermite interpolation properties between the maximum singular values functions of the full and the reduced transfer functions.This problem can be solved by evaluating the full and the reduced transfer function at the interpolation points, checking whether the Hermite interpolation properties are satisfied, and expand the projection spaces by including further singular vectors if necessary.Remarkably, in practice it happens very rarely that the Hermite interpolation conditions are not satisfied.
Also note that in course of the iteration we use i ω j as a further initial interpolation point for warmstarting the computation oof F (g j+1 , •) H∞ in line 8 of Algorithm 1 for faster convergence.Moreover, the value of g (j) is used as an initial value for the optimization of g → F j (g, •) L∞ in line 7 of Algorithm 1.

Choice of the Initial Interpolation Points
Recall that our method for computing the optimal gains g * is based on Hermite interpolation between the original and the reduced parameter-dependent transfer functions.Thus, in case of convergence, we would only know that ( 6) is satisfied locally.While convergence to a local minimizer g * of g → F (g, •) H∞ would be feasible, convergence to a local maximizer of ω → F (g, iω) 2 would be troublesome, since this would underestimate the closed-loop performance for a given g.
For this reason it is important to find good initial points in line 2 of Algorithm (1) to enhance the chance of converging to a global maximizer of the singular value function.Thus, for g (i) , . . ., g (ν) , we estimate the location of the dominant poles of the transfer function which are responsible for large maximum singular values of the transfer functions F ( g (i) , •) on the imaginary axis.
According to [27,5], a transfer function m× can be written in pole/residue representation as where λ i and R i , i = 1, . . ., 2n, are the poles and the residues of the transfer function, respectively.Now let (λ Definition 5.For the transfer function The k most dominant poles are the k poles with the largest values of The dominant poles have a close relationship to the local maxima of the function ω → F (iω) 2 .Namely, if λ i is a dominant pole of F (s), then it is likely that a large local maximum can be found near Im(λ i ).Conversely, if ω → F (iω) 2 has a large local maximum at ω 0 , then a dominant pole will likely be close to iω 0 .A graphical interpretation of this relationship is available in [20, p. 17].A similar notion of pole dominance has already been used in the context of damping optimization, see [5].
In principal, one could compute a number of dominant poles of F ( g (i) , s) ∈ R(s) m× for all g (i) , i = 1, . . ., ν and use their imaginary parts as initial frequencies ω i,j in line 2 of Algorithm 1.This can be done using the subspace accelerated MIMO dominant pole algorithm (SAMDP) [27].
On the other hand, if the state-space dimension n is of moderate size, then we can also estimate R i 2 and λ i by transforming the system to modal coordinates before.This will be outlined next.We estimate R i 2 using the transfer function of the undamped problem which is Since M and K are positive definite matrices, there exists a matrix Φ = φ 1 . . .φ n ∈ R n×n such that Therefore, we can write Then by (14), for i = 1, . . ., n, the residues of F 0 are given by where e i ∈ R n is the ith unit vector in R n and γ i is a scaling parameter such that This gives and hence, we get Next we estimate the pole locations of the transfer function F (s).By a linearization of the quadratic eigenvalue problem the poles of F (s) are eigenvalues of the matrix Here we consider the matrix A as a perturbation of the matrix Note that the eigenvalues of A 0 are λ ± 0,i = ±iω i with associated right and left normalized eigenvectors (both are the same) .
By a standard result from first order perturbation theory [28], it holds that the eigenvalues λ ± i , i = 1, . . ., n of A are given by which gives the approximation This approximation is usually good, if the norm of the perturbation given by C(g) 2 is small.Finally we sort ω i according to the values of in decreasing order, which can be done very cheaply once the system is known in modal coordinates.

Numerical Experiments
In this section we consider numerical experiments.In order to illustrate the efficiency of our proposed approach, we compare the optimal gains calculated for the full-order system with the optimal gains obtained by Algorithm 1.In each of our experiments we use the open-source package GRANSO for the parameter optimization.Alternatively, one could use HANSO § , see also [14,22].However, since we wish the optimal parameters to be in a feasible domain, we have a constrained optimization problem which makes the use GRANSO more appropriate in our situation.For the full-order problem, the H ∞ -norms are calculated by the Boyd-Balakrishnan algorithm for descriptor systems [6] which is implemented in the FORTRAN routine AB13HD.f and is called by the mex file linorm h.f in MATLAB.On the other hand, in Algorithm 1, we compute the H ∞ -norms for a sequence of large-scale non-parametric problems by employing linorm subsp and then optimize the gains on a sequence of reduced parametric systems.
In our experiments, the computations have been carried out on a machine with four Intel R Core TM i5-4590 CPUs @ 3.30 GHz and 16 GB RAM.The results reported in this work are calculated by MATLAB version 9.6.0.1072779 (R2019a) on a 64-bit Linux operating system.
Example 1.We consider an n-mass oscillator which describes the mechanical system of n masses and n + 1 springs shown in Figure 1.The mathematical model is given by (1a)-(1c), while the mass and stiffness matrices are The masses and stiffnesses have the following configuration: We are interested in the states that correspond to the masses with indices ranging from 290 to 309, that is, in this example we consider the performance output . . . .
, where the internal damping C int is given by (2) and where we choose the critical damping model which is widely used in the literature, see, e. g., [12,8,5].Note that this choice makes the unforced system (1a) asymptotically stable.
We consider two dampers with different gains, that is, we choose the matrix G = diag(g 1 , g 2 ) ∈ R 2×2 .The geometry of the external damping depends on the dampers' positions (j, k) which are encoded in the matrix B 2 by setting B 2 = e j − e j+1 e k − e k+1 , where e j and e k are the jth and the kth canonical vector in R n .
In general, the problem is to optimize the positions of the dampers, but this requires an optimization of the gains for a large number of different damper positions.Thus, we illustrate the efficiency of our method for the optimization of the damper positions (j, k) for j ∈ {40, 140, 240, 340, 440, 540} and k ∈ {60, 160, 260, 360, 460, 560}, so in total we obtain 36 different configurations for which we have to perform a gain optimization.
Moreover, the parameter α c determines the influence of an internal damping as shown in (2).The internal damping has a strong influence on the system, so we show the results of our approach for two different settings: Problem a) α c = 10 −5 : In this case we have a very small internal damping which means that almost all eigenvalues of the matrix polynomial s 2 M + sC(0) + K (that corresponds to the problem without dampers) are very close to the imaginary axis.The optimization of this problem is particularly difficult, since these eigenvalues often introduce a very large number of thin peaks (local maxima) in the function ω → σ max (F (g, iω)).
Problem b) α c = 10 −2 : In this case the internal damping is moderate which improves asymptotic stability.Moreover, the amount of large local maxima in ω → σ max (F (g, iω)) is moderate which makes it easier to compute the H ∞ -norm in this case.This choice of α c is more realistic from the practitioner's point of view.
For both cases we use the following general computational setup and parameters in Algorithm 1: • The tolerance of GRANSO for (approximate) stationarity is set to 10 −12 .
• We use linorm subsp with default parameters and with tangential interpolation.For each call of linorm subsp, the 30 most dominant poles computed according to Subsection 3.2 are chosen as initial interpolation points.In case we use SAMDP, we first compute 180 dominant poles and pick the 30 most dominant ones out of these.
• Initially, we choose 4 initial parameters to set up the initial reduced-order model as follows.
• The relative termination tolerances for the gains and the L ∞ -norm are both set to 10 −6 .In other words, we terminate, if

2
, or or the maximum number of iterations (which we have set to 30) has been reached.
We further test our method with different computational modes whose results we present further below: Mode i) For each g (1) , . . ., g (4) , we perform tangential interpolation of F ( g (i) , •) at 30 equidistantly distributed samples i ω i,1 , . . ., i ω i,30 in the frequency range ω i,j ∈ [0, ω max ] for j = 1, . . ., 30 to setup the initial reduced-order model.Here ω max is the maximum modulus among all eigenvalues of the matrix polynomial s 2 M + K. To determine initial interpolation points for linorm subsp we choose the heuristic approach from Section 3.2.
Mode ii) We construct the initial reduced-order model as in Mode i).To determine initial interpolation points for linorm subsp we use the SAMPD algorithm.
Mode iv) We construct the initial reduced-order model as in Mode iii).To determine initial interpolation points for linorm subsp we use the SAMDP algorithm.
In the following we compare our results to the naive optimization approach that consists of optimizing the full-order problem with GRANSO and the initial point set to be the first initial parameter, that is the initial point for Problem a) is ( 1010 ), while for Problem b) it is set to ( 100 100 ).Figures 2 and 4 show the relative errors in the computed optimal gains for Problem a) and Problem b), respectively, where as reference values we choose the optimal gains computed for the full-order model.More precisely, for the ith configuration, the Figures display the values of g Since the objective function may be flat near the optimizer, i. e., there may be a large area of almost optimal gains, the differences in the optimal values of the objective function are more informative when assessing the quality of the results.
As expected, our method has some difficulties with Problem a), since for a very small internal damping, the mappings ω → F (g, iω) 2 have a huge number of local maxima and then it may e. g. happen that the global maximum is missed by our method.Further, these local maxima often lead to extremely steep peaks in the maximum singular value plots which can cause additional numerical difficulties since then the maximum singular values of the transfer function are very sensitive with respect to small changes in ω.These problems can be seen in Figure 3, where some of the errors are quite large.This behavior is a clear limitation of Algorithm 1.The difficulty of this problem is further illustrated in Figure 6 in which we plot F (g * , iω) 2 as well as F r g * , iω 2 over ω (where F r denotes the final reduced transfer function and g * are the optimal gains computed by Algorithm 1.) for one of the configurations of Problem a).Remarkably, even if the original transfer function induces extremely many peaks in its maximum singular value plot, this is not the case for the reduced one, but the two global maximizers still coincide approximately.
On the other hand, as shown in Figure 5, the results are more convincing for a moderate internal damping which is also more realistic from the application's point of view.
We also compare the computational time needed for the optimization process.For each configuration we measure the runtimes needed to optimize the full-order model with the naive approach and the   1.From this table we can conclude that Algorithm 1 accelerates the optimization process considerably.For our examples, the initial interpolation points can be calculated faster by using formula (17) compared to using the SAMDP algorithm, i. e., see the entries corresponding to Mode i) and Mode iii) in Table 1.However, formula (17) can no longer be efficiently applied, if the system dimension becomes very large.On the other hand, the SAMDP algorithm could still potentially be used in this case, since only a small fraction of the eigenvalues may be computed.For this, it is necessary to be able to solve the involved linear systems efficiently.For instance, this is the case if instead of critical damping one uses Rayleigh damping, since then C int = αM + βK for some small constants α, β > 0 would still be a sparse matrix.We do not consider this case here, since for very large systems, the computations for the naive optimization approach become prohibitively expensive.From the figures we can further see that it is more efficient to evaluate the H ∞ -norms at the initial gains in order to set up the initial reduced order model (Mode iii) and Mode iv)) compared the sampling procedure used in Mode i) and Mode ii).
We have also evaluated our approach with full instead of tangential interpolation.In our experience, this does not have a significant influence on the number of iterations but evaluating the H ∞ -norm of the reduced transfer functions becomes very expensive, since the projection space dimensions grow very fast.Therefore, in this case the runtime is in the same order of magnitude as for the naive method.

Conclusion
Altogether we can conclude that with our new approach we have been able to perform the semi-active H ∞ damping optimization for a problem with moderate internal damping with satisfactory relative accuracy, while the optimization process was considerably accelerated.On the other hand, our method has a few problems with some of the configurations with very small internal damping.Such problems must be treated more carefully and it is necessary to have a mindful choice of the initial sampling data as well as of the algorithm parameters.Still, the optimization problem could be solved with a satisfactory accuracy for most of the configurations of this extremely hard problem.

Figure 1 :
Figure 1: An n-mass oscillator with two dampers * and g (i) * denote the gains obtained by optimizing the full-order model directly and by applying Algorithm 1, respectively.Figures 3 and 5 on the other hand show the relative errors in the computed H ∞ -norms.In other words, the Figures present the values of | F (g (i) * ) H∞ − F ( g (i) * ) H∞ |/ F (g (i) * ) H∞ for the ith configuration.

Figure 2 :
Figure 2: Comparison of optimal gains for Problem a)

Figure 3 :
Figure 3: Comparison of the optimal H ∞ -norm for Problem a)

Table 1 :
Average time ratio