## 4 Identification Procedure  
The term _pole-zero identification_ refers to the achievement of the poles and zeros of a linear, or linearized, system described by its frequency response. This is usually done using optimization techniques able to fit a given frequency response of the linear system towards a transfer function defined as the ratio of two polynomials [36; 37]. This kind of linear system identification in the frequency domain has several applications in a wide variety of engineering fields and it can be constrained or unconstrained. In this work, the pole-zero identification is performed adopting a constrained procedure taking advantage from the \"fminmax\" Matlab function that seeks the minimum of a nonlinear error function \(W\) of n real variables, the Zeros and the Poles. For the present identification problem, this function will be defined as the weighted sum of the normalized root mean square error (NRMSE) calculated for the storage modulus and the loss tangent as in (11):  
\[W=(1-r)\cdot Err_{E}+r\cdot Err_{tan5} \tag{11}\]  
where r is the weight factor and \(Err_{tan5}\) and \(Err_{E}\) are defined by (12)  
\[Err_{tan5}=\frac{\sqrt{\frac{\sum_{i=1}^{N}\left(tan5-tan5_{model}\right)^{2}} {N}}}{mean(tan5)}\quad Err_{E}=\frac{\sqrt{\frac{\sum_{i=1}^{N}\left(E-E_{ model}\right)^{2}}{N}}}{mean(E)} \tag{12}\]  
Notice that this way of defining the error function allows to fit both modulus and phase at the same time. As concerns the weighting coefficient \(r\), generally, it can be set equal to the more convenient value in order to have the more robust results. In the present paper, the weighting coefficient \(r\) has been assumed equal to 0.5 for all identifications in order to compare the results of different procedures. The function \(W\) must be constrained imposing that the poles \(p\) and zeros \(z\) are alternate, according to the third property demonstrated by [PERSON] [38] for linear dissipative systems (13):  
\[z_{1}<p_{1}<z_{2}<p_{2}<z_{3}<p_{3}\ldots<z_{i}<p_{i} \tag{13}\]  
As concerns the initial condition of the identification procedure, the initialization of the pole-zero parameters has been carried out by means the method proposed by [PERSON] et al. [32]. Dividing the frequency domain into N equal frequency subdomains each Pole Zero couple is given by (14):  
\[\left\{ \begin{array}{l}\omega_{z,i}=(1-r)\cdot\omega_{z,i}^{Modulus}+r\cdot \omega_{z,i}^{Phase}\\ \omega_{p,i}=(1-r)\cdot\omega_{p,i}^{Modulus}+r\cdot\omega_{p,i}^{Phase}\end{array} \right. \tag{14}\]where, posing \(\chi=\log_{10}(\omega)\), \(\omega_{z,i}\) and \(\omega_{p,i}\) for modulus and phase are calculated as in (15):  
\[\left\{ \begin{array}{l}\chi_{z,i}^{Modulus}=\chi_{c,i}-\frac{\alpha_{i}}{2}\\ \chi_{P,i}^{Modulus}=\chi_{c,i}+\frac{\alpha_{i}}{2}\end{array} \right.\ ;\ \ \ \left\{ \begin{array}{l}\chi_{z,i}^{Phase}=\chi_{c,i}-\frac{\alpha_{i}}{2}\\ \chi_{P,i}^{Phase}=\chi_{c,i}+\frac{\alpha_{i}}{2}\end{array} \right. \tag{15}\]  
where \(\chi_{c,i}\) is the medium frequency of a couple and \(\alpha_{i}\) is evaluated considering the stiffing of the modulus \(|E^{*}(\omega)|\) and the area under the phase curve the \(\delta(\omega)\) defined respectively as in (16) and (17):  
\[S=\alpha_{i}=\int_{\chi_{q},i}^{\chi_{0,i}}\frac{d}{d\chi}\left(\log_{10}\big{|} E*(\chi)\big{|}\right)d\chi \tag{16}\]  
\[A=\alpha_{i}\frac{\pi}{2}=\int_{\chi_{q},i}^{\chi_{0,i}}\delta(\chi)d\chi \tag{17}\]  
From the computational point of view, the convenience of the pole zero formulation in the identification of parameters of transfer functions lies in how the boundary conditions and the initial conditions of the optimization problem can be defined. As an example, using the formulation (6) of the GM model, the parameters to be identified are stiffness and damping coefficients. Hence, defining both the initial conditions and the boundary conditions is a hard task because they can vary between zero and very large numbers and as a function of the elements of the model itself.