## 2 Implementation of the Nevanlinna-Pick Constraint into Structural Dynamics  
According to the boundary NP interpolation theorem [48], the existence of a solution, which is a rational function mapping from the open right half plane \(\pi^{+}\) into an open unit disk \(D\), is related to the positive definiteness of the following Pick matrix \(\mathbf{\Lambda}_{lm}\)  
\[\Lambda_{lm}=\left\{ \begin{array}{c}\frac{1-\overline{f}_{l}f_{m}}{\overline {s}_{l}+s_{m}}\quad\quad,\quad if\quad l\
eq m\\ \left|f^{\prime}(s_{l})\right|\quad,\quad if\quad l=m\end{array} \right. \tag{1}\]  
where \(f_{k}\) is the given complex function defined by  
\[\left|f_{k}\right|=1 \tag{2}\]  
and \(s_{k}=i\omega_{k}\) is the complex Laplace transform variable, \(i\) is the imaginary unit, and \(\omega_{k}\) is the excitation frequency. In summary, a physical system cannot display an arbitrary frequency response function (FRF) \(f(i\omega)\) and derivative \(f^{\prime}(i\omega)\) at discrete frequencies \(\omega_{k}\). The above equations are implemented into finite element equations for sample structural problems, represented by the following equation:  
\[M\ddot{x}+C\dot{x}+Kx=F\]where \(M\), \(C\), and \(K\) are, respectively, mass, damping, and stiffness matrices. \(F\) is the applied force vector, \(x\) is the displacement vector, and time rates are denoted by superposed dots. If the Laplace transform of Equation (3) is taken, then we obtain the following equation in terms of the Laplace transformed variables:  
\[\Big{(}M\!s^{2}+Cs+K\Big{)}X(s)=F(s) \tag{4}\]  
where \(s\) is the Laplace transform variable. The corresponding system transfer function is expressed as follows:  
\[G(s)=\frac{F(s)}{X(s)}=\Big{(}M\!s^{2}+Cs+K\Big{)}^{-1} \tag{5}\]  
and \(f_{k}=G(\omega_{k})\). For a harmonic force with the frequency \(\omega\), i.e., \(F(t)=F_{0}e^{i\omega t}\), the steady-state response of a structure can be expressed as follows [1]:  
\[x_{p}(t)=x_{0}|G(i\omega)|e^{i(\omega t+\varnothing)} \tag{6}\]  
where \(x_{p}(t)\) denotes the particular solution, \(|G(i\omega)|\) is the magnitude of the transfer function of the structure, and \(\varnothing\) is the phase angle. The NP equations are thus related to the steady-state response of the structure to a harmonic force. To derive insight into imposing the NP constraint equations onto structural equations, the physical meaning of the terms in Equation (1) is pursued. To this end, following Equation (5), it can be written that \(X(i\omega)=G(i\omega)F(i\omega)\). Suppose that the forcing term in Equation (3) is an impulsive force with the following harmonic expansion: \(F(t)=F_{0}\delta(t)=F_{0}\sum_{i=1}^{\infty}A_{i}\text{sin}\left(\omega_{i}t\right)\), where \(F_{0}\) is a constant vector, \(\delta(t)\) is the impulsive force, and \(\omega_{i}\) the are harmonic frequencies. Then, the impulse response \(g(t)\) can be calculated by using the Bode gain and phase equations [1]:  
\[g(t)=\sum_{i=1}^{\infty}F_{0}|G(i\omega_{i})|A_{i}\text{sin}(\omega_{i}t+ \varphi_{i}) \tag{7}\]  
where \(\varphi_{i}=Arg(G(\omega_{i}))\) is the phase angle and \(G(i\omega)\) is found using the Laplace transform of \(g(t)\) as follows: \(G(i\omega)=\mathcal{L}\{g(t)\}\). On the other hand, since \(G^{\prime}(i\omega)=-\mathcal{L}\{tg(t)\}\) and \(|G^{\prime}(i\omega)|=|\mathcal{L}\{tg(t)\}|\), e.g., see [1], we derive  
\[tg(t)=\sum_{i=1}^{\infty}F_{0}|G^{\prime}(\omega_{i})\big{|}A_{i}\text{sin}( \omega_{i}t+\varphi_{i}{}^{*}) \tag{8}\]  
where \(\varphi_{i}{}^{*}=Arg(G^{\prime}(\omega_{i}))\) are phase angles. In summary, Equations (7) and (8) indicate that the frequency domain quantities of \(|G(\omega_{i})|\) and \(|G^{\prime}(\omega_{i})|\) are related to the impulse response magnitude and time-weighted impulse response magnitude, respectively. As a consequence, the magnitudes of \(g(t)\) and \(tg(t)\) cannot be independent due to Equation (1). Consequently, imposing the NP constraint Equation (1) in optimization costs assumes the dependence between \(g(t)\) and \(tg(t)\).