© The Author(s), Published by EDP Sciences, 2025

1 Motivation

The application of field simulation tools, such as the finite element method, is now commonplace in many design processes. The modelled, virtual prototypes are used to optimise design goals before physical components are manufactured to reduce resource usage and development time. A prerequisite for a simulation-driven design is that the simulation accurately represents the physical behaviour of the component, thus requiring exact quantitative knowledge of the underlying material parameters. These parameters are especially hard to determine if the physical process relies on the interaction of different phenomena, for example in the piezoelectric effect, where the mechanical, electrical, and coupling behaviour needs to be quantified. Design processes of piezoelectric devices, such as ultrasound transducers, wire bonding machines, and dental scalers thus require of full, quantitative material description to accurately describe the electromechanical conversion processes thatoccur.

A standardised method to determine a full set of piezoceramic material parameters is given in the IEEE Standard on Piezoelectricity [1], which proposes measurements of the resonance frequencies of four different samples (a long cylinder, a flat disc, and two differently polarised cuboids) of the same material. The standard provides analytical expressions that can be used to determine the material parameters from the resonance and anti-resonance frequencies of each sample. This procedure can, however, lead to inconsistent material parameter vectors, due to inherently different processing steps for each of the four samples. Further research on the determination of the material parameters of piezoceramics acknowledges the fact that analytical relations yield only approximative results: Recent approaches for quantitative material modelling primarily rely on solving inverse problems and matching the physical behaviour of samples with simulations [2, 3]. A prerequisite for the application of these methods is that the quantity to be matched by the simulation, usually the frequency-dependent, electrical impedance of the samples, is sufficiently sensitive with respect to all material parameters. Kybartas et al. [4] show, that this can be achieved by choosing a sample geometry that promotes coupled resonance modes. Alternatively, Rupitsch et al. [5] propose an inverse method for a complete piezoceramic characterisation that employs two cuboid samples and incorporates measurements of the surface velocity alongside the electrical impedance measurement. Aiming to perform the material characterisation using a sample geometry identical to the target application, efforts are made to optimise the sensitivity to all material parameters by structuring the electrodes of disc-shaped samples [6–9]. These modifications allow the determination of a full set of piezoelectric material parameters, including different damping models, using only a single sample [10].

In high-power acoustic applications, piezoceramic elements with a ring geometry are usually applied. For practical applications, the sensitivity of the impedance with respect to the material parameters is considered sufficient to determine a full set of piezoelectric material parameters using an unmodified ring-shaped specimen [11].

The solution procedure for the inverse problem arising when trying to match the measured electrical impedance with the simulation model output requires multi-stage approaches [8, 10, 11]. Initial values are estimated by evaluating approximative analytical relations derived from the IEEE Standard on Piezoelectricity [1] and complex expressions for the electrical impedance [12]. The results of the estimation are subjected to an intermediate optimisation step, modifying only single parameters based on particularly sensitive resonance frequencies. The main optimisation comprises a complex, sensitivity-based, block coordinate descent procedure that needs to be adapted to the sample geometry at hand [10].

The primary aim of this study is to create a simplified, more robust material parameter identification method using data-driven approaches e.g. for an improved initial value estimation method. Results not only show a more straightforward optimisation strategy, and thus reduced reliance on empirical approaches in designing said strategy, but also a significantly reduced computational effort for the identification procedure.

The structure of this article is as follows: In Section 2 the material description and the finite element (FE) simulation model used in the identification procedure are presented. This model is used to generate synthetic training data (Sect. 3). An approach to detect artefacts occurring in physical measurement systems based on the training data is also described. The multi-stage inverse procedure for material parameter identification is described in Section 4, including initial estimates generated by a neural network trained as the approximate inverse of the FE model, and a fast, intermediate optimisation step using a neural network acting as a surrogate of the model. A final optimisation step based on the FE model yields results for the material parameters, which are presented in Section 5 for different sample temperatures.

2 Numerical model

To describe the linear behaviour of a piezoelectric material, mechanical and electrical properties as well as their coupling need to be quantified. For the mechanical part Hooke’s law describes the relation between stress T and strain S inside the material via the elasticity matrix c^E. This description is expanded to include the effects of piezoelectricity and thus the influence of the electric field E [13]:

$$\begin{array}{c}\hfill \mathit{T}={\mathit{c}}^{\mathrm{E}}\mathit{S}-{\mathit{e}}^{\mathrm{t}}\mathit{E},\end{array}$$(1)

where e^t is the transposed piezoelectric coupling matrix. The defining equation for the dielectric displacement D is similarly supplemented by a term that describes the coupling to the mechanical strain S [13]:

$$\begin{array}{c}\hfill \mathit{D}={\mathit{\epsilon }}^{\mathrm{S}}\mathit{E}+\mathit{e}\mathit{S},\end{array}$$(2)

with ε^S denoting the electric permittivity at constant strain.

The samples examined in this study are polarised, piezoelectric ceramics, which can be considered transversely isotropic [1]. With the polarisation direction in 3-direction, the material matrices in Voigt notation [14] are occupied as follows:

$$\begin{array}{c}\hfill {\mathit{c}}^{\mathrm{E}}=\left[\begin{array}{cccccc}{c}_{11}^{\mathrm{E}}& c_{12}^{\mathrm{E}}& c_{13}^{\mathrm{E}}& 0& 0& 0\\ [3pt]c_{12}^{\mathrm{E}}& c_{11}^{\mathrm{E}}& c_{13}^{\mathrm{E}}& 0& 0& 0\\ [3pt]c_{13}^{\mathrm{E}}& c_{13}^{\mathrm{E}}& c_{33}^{\mathrm{E}}& 0& 0& 0\\ [3pt]0& 0& 0& c_{44}^{\mathrm{E}}& 0& 0\\ [3pt]0& 0& 0& 0& c_{44}^{\mathrm{E}}& 0\\ [3pt]0& 0& 0& 0& 0& \frac{c_{11}^{\mathrm{E}}-c_{12}^{\mathrm{E}}}{2}\end{array}\right]\end{array}$$(3)

for the elasticity matrix,

$$\begin{array}{c}\hfill {\mathit{\epsilon }}^{\mathrm{S}}=\left[\begin{array}{ccc}{\epsilon }_{11}^{\mathrm{S}}& 0& 0\\ [3pt]0& {\epsilon }_{11}^{\mathrm{S}}& 0\\ [3pt]0& 0& {\epsilon }_{33}^{\mathrm{S}}\end{array}\right]\end{array}$$(4)

for electric permittivity matrix, and

$$\begin{array}{c}\hfill \mathit{e}=\left[\begin{array}{cccccc}0& 0& 0& 0& e_{15}& 0\\ [3pt]0& 0& 0& e_{15}& 0& 0\\ [3pt]e_{31}& e_{31}& e_{33}& 0& 0& 0\end{array}\right]\end{array}$$(5)

for the piezoelectric coupling matrix, yielding ten independent material parameters [13]. For the present study, mechanical losses are quantified by including a Rayleigh damping model [15], adding two degrees of freedom: Damping proportional to mass (α_M), which models absorption decreasing with frequency, and damping proportional to stiffness (α_K), which models absorption increasing with frequency. Electrical and coupling losses are assumed to be negligible. For a piezoelectric material with the description assumed here, a parameter vector p with twelve elements must be identified:

$$\begin{array}{cc}\hfill \mathit{p}=& [c_{11}^{\mathrm{E}},c_{12}^{\mathrm{E}},c_{13}^{\mathrm{E}},c_{33}^{\mathrm{E}},c_{44}^{\mathrm{E}},{\epsilon }_{11}^{\mathrm{S}},{\epsilon }_{33}^{\mathrm{S}},\hfill \\ \hfill & e_{15},e_{31},e_{33},{\alpha }_{\mathrm{M}},{\alpha }_{\mathrm{K}}].\hfill \end{array}$$(6)

The density of the sample ρ is directly accessible through geometric and gravimetric measurements, and thus not an element of the target parameter vector.

Using this material description, a numerical finite element model is implemented in openCFS [16] with geometry matching the physical samples analysed in this study: a piezoelectric ring with a thickness t of 1 mm, an outer radius r_o of 6.35 mm, and an inner radius r_i of 2.6 mm (Fig. 1). The model is publicly available as part of the training data set [17]. Due to the rotational symmetry, the simulation model can be simplified to the cross-sectional area, significantly reducing the computational effort. This area is meshed with an irregular mesh with an approximate element size of 200 μm and third order polynomials. The elastic wave equation in frequency domain and the electrostatic equation are solved for the cross-section. The electrodes present in the physical sample are realised as constrains for the electric field at the upper and lower face of the ring. Oscillations are induced by an electric potential U_sim applied through boundary conditions of the electric field at the upped face of the ring, while the potential at the bottom face is fixed to zero. The typical output of this model is the frequency domain solution for the mechanical displacement and the electric potential, however, of primary interest for this study is the resulting charge on the electrodes Q_sim. Using the definition for the electric impedance Z_sim and the relationship between electric current I_sim and charge in thefrequency domain, the model can be used to determine the frequency-dependent impedance of the virtual sample:

Figure 1 Sectional view of the simulation model of a piezoelectric ring used in the present study. Due to rotational symmetry, a two-dimensional mesh of the cross-section area A can be used.

$$\begin{array}{c}\hfill Z_{\mathrm{sim}}=\frac{U_{\mathrm{sim}}}{I_{\mathrm{sim}}}=\frac{U_{\mathrm{sim}}}{\mathrm{j}\omega Q_{\mathrm{sim}}},\end{array}$$(7)

where j is the imaginary unit and ω is the angular frequency. The model thus provides the frequency-dependent electrical impedance vector Z_sim for a given set of material parameters p. Examples for the electrical impedance data are shown in Section 3.1, with local minima and maxima indicating resonance and anti-resonance frequencies, respectively [13]. However, the task of determining the material parameters for a given electrical impedance is not trivial and central to this study.

2.1 Sensitivity analysis

Several studies explore the sensitivity of the electrical impedance of piezoelectric samples with respect to the material parameters [5–7, 11, 18]. Analyses are conducted by varying each material parameter p_i by ±1% in the simulation model and observing changes in the squared logarithmic absolute value of the electrical impedance Z_sim:

$$\begin{array}{c}\hfill {{\rm Y}}_i=\sum {\left({log}_{10}|{\mathit{Z}}_{\mathbf{sim},p_i+1\%}|-{log}_{10}|{\mathit{Z}}_{\mathbf{sim},p_i-1\%}|\right)}^2.\end{array}$$(8)

An initial data set for the material parameter is required for the sensitivity analysis, which is assumed as: $c_{11}^{\mathrm{E}}=130\mathrm{GPa}$, $c_{12}^{\mathrm{E}}=84\mathrm{GPa}$, $c_{13}^{\mathrm{E}}=83\mathrm{GPa}$, $c_{33}^{\mathrm{E}}=117\mathrm{GPa}$, ${\epsilon }_{11}^{\mathrm{S}}=8.15{\mathrm{nF}\mathrm{m}}^{-1}$, ${\epsilon }_{33}^{\mathrm{S}}=6.58{\mathrm{nF}\mathrm{m}}^{-1}$, $e_{15}=12.1{\mathrm{C}\mathrm{m}}^{-2}$, $e_{31}=-6.03{\mathrm{C}\mathrm{m}}^{-2}$, $e_{33}=15.5{\mathrm{C}\mathrm{m}}^{-2}$, ${\alpha }_{\mathrm{M}}=120{\mathrm{ms}}^{-1}$, and ${\alpha }_{\mathrm{K}}=0.1\mathrm{ns}$ with a density ρ of 7800 kg m⁻³. Results for the sensitivity of the electrical impedance with respect to changes in the material parameters are visualised in Figure 2 for piezoelectric discs with electrodes fully covering the faces, the piezoelectric rings central to this study, and for piezoelectric discs with electrodes structured for optimised sensitivity [8]. The parameters are orderer by the results for the piezoelectric disc, which shows that the sensitivity of the electrical impedance with respect to some material parameters is particularly small for quantities describing the lateral behaviour like ${\epsilon }_{11}^{\mathrm{S}}$ or e₁₅. To remedy this, methods to determine an optimised electrode structure to increase said sensitivity are developed [7–9], the result of which are clearly visible in the sensitivity analysis (Fig. 2). The sensitivities of the piezoelectric ring geometry considered in this study fall in-between the results for the other sample geometries considered. While the sensitivities are certainly less advantageous compared to a geometry optimised for sensitivity [8], they are considered adequate for a material characterisation [11]. Additional optimisation of the piezoelectric ring’s electrode structure to increase sensitivity is thus not considered necessary.

Figure 2 Results of a sensitivity analysis of the numerical model for piezoelectric discs with full electrodes (Full), the piezoelectric rings (Rings) studied here [17], and the piezoelectric discs with structured electrodes optimised for sensitivity (Structured) [8, 19]. Material parameters are varied by ±1% and the squared change in logarithmic absolute value of the electrical impedance is observed.

Overall, the analysis indicates large differences in sensitivity for the different parameters, covering several orders of magnitude. This can be disadvantageous for the subsequent optimisation procedure and is addressed in several studies, as previously mentioned [7–9, 11, 18]. However, the numerical values shown here hide the fact that the influence of each parameter on the electrical impedance is qualitatively different. The elastic parameter $c_{12}^{\mathrm{E}}$, for example, is known to only affect the first radial resonance frequency in axisymmetric samples [18]. The damping parameter α_K, which show the lowest sensitivity for all considered sample geometries, influences the damping of all higher frequency resonance modes and can typically be determined reliably, even if the sensitivity is low [10]. In turn, values of the damping parameters α_M and α_K for piezoceramics cover wider ranges than the other material parameters.

The qualitative difference in the influence that each parameter has on the electric impedance of the virtual sample can also be used to diagnose possible misidentified parameters [20]. This is especially relevant because it is as of now unknown if the objective function of piezoelectric material parameter identification problem based on the electrical impedance has a single, prominent, local minimum. A recent study did, however, show that different initial value combinations yield very similar optimisation results [21]. Unsuccessful optimisation processes are typically easily identifiable by high objective function values and a clear mismatch between virtual Z_sim and measured impedance Z_meas, with the character of themismatch indicating which of the parameters are incorrectly identified. Cases in which different material parameter value combinations yield a very similar virtual electrical impedance Z_sim, which would indicate multiple solutions to the identification problem, have, to the authors’ knowledge, not been reported.

3 Training data generation and analysis

A prerequisite for the application of machine learning methods is the availability of sufficient training data. For the specific application discussed here, this amounts to numerous data sets composed of the sampled, complex, frequency-dependent impedance of a piezoelectric ceramic and the associated material parameters. For this study, the geometry of the piezoelectric ceramic is considered constant. Data sets of this configuration that are determined from physical measurements are currently not available in large numbers, necessitating the generation of synthetic training data. This data is generated by evaluating a finite element simulation model [16] of the sample with randomised material parameter. The model is also applied in the subsequent optimisation procedure to solve the inverse problem for material parameter determination and is thus already designed for short evaluation times [10]. The utilisation of synthetic training data implies that the physical behaviour of the sample can be reproduced by the simulation model. This assumption is, however, inherent to the inverse measurementprocedure.

3.1 Training data generation

To generate synthetic training data, the material parameter values are independently sampled from continuous uniform distributions with intervals [a, b] (Tab. 1). The intervals are chosen based on results from literature [5, 10, 18] with some margin added. The damping parameters are sampled from logarithmic distributions to account for the larger interval. The simulations are performed in frequency domain with 1999 samples up to 8 MHz and yield the complex electrical impedance vector Z_sim. It is assumed that at least 2¹² to 3¹² data sets are required to train a network for the estimation of 12 parameters. If sampling the parameter space with a regular grid, this would amount to 2 to 3 samples in each direction respectively. In terms of polynomials, this would only allow for the identification of linear or quadratic expressions for each parameter and is thus still considered a sparse sampling of the parameter space.

As a trade-off between the considered bounds for the number of samples and economical use of computing resources, 300 000 data sets are generated by evaluating the model in parallel on a local computational cluster, yielding an equal number of synthetic electrical impedances with the corresponding piezoelectric material data. Four randomly selected examples for the impedance data are visualised in Figure 3, showing characteristic features such as an overall capacitive behaviour and pairs of resonance and anti-resonance frequencies indicated by local minima and maxima respectively. The results are compiled in a database and are publiclyaccessible [17].

Figure 3 Randomly selected examples for the synthetic impedance data from the data set [17].

3.2 Measurement artefact detection

The large database can also be utilised to determine if a given measurement data set is within the scope of the training data or to detect if the measurement may be influenced by artefacts. To realise this detection method, a simple, fully connected autoencoder network [22] with twelve latent variables (equal to the number of parameters) is trained. The setup of the autoencoder is symmetrically with input and output layers of 2000 neurons (1999 impedance samples and the density ρ) and four hidden layers in the encoder (1000 neurons, 500 neurons, 300 neurons, and 100 neurons) and decoder. Hyperbolic tangent activation functions are used for all hidden layers. The output and the latent space layer are set up with exponential linear units to enable larger numerical values. The logarithm of the absolute values of the impedance is supplied to the autoencoder during training and the density is scaled in 1 g cm⁻³ = 1000 kg m⁻³. Adaptive momentum estimation [23] is used to train the autoencoder for 1000 epochs. For the subsequent analysis, the numerical values of the latent variables are observed when supplying the training data set to the autoencoder. Even though the number of neurons in the hidden layer is equal to the number of parameters, the values of the latent variables can not be interpreted physically but are analysed for their statistical properties.

It is assumed that measurements, which can be modelled by determining adequate material parameters, yield values of the latent variables that are within the distribution yielded by the training data set. As an example, physical measurement data from piezoelectric rings made of two different materials (PIC181 and PIC184, PI Ceramic) are analysed (Fig. 4). The box plot shows the distributions resulting from the training data set and the markers indicate the values yielded by the encoder when applying measurement data. It is observable that the values of the latent variables of both measurement data sets are within the range of the distributions of the training data.

Figure 4 Distribution of the values in the latent variables of the trained autoencoder for the training data set (box plot) and values for specific measurements of different piezoelectric ceramics (PIC181 and PIC184, PI Ceramic) plain and with superimposed artificial measurement artefacts.

As examples, synthetic measurement artefacts are superimposed to the impedances in the form of additional capacitive and inductive impedances. These may occur in physical measurement setups due to imperfections in contacting the sample or long electrical lines. It can be observed (Fig. 4) that the resulting values of the latent variables are at least partially found outside the distributions yielded from the training data set. This deviation from the expected result for unaffected measurements can be expressed numerically. The average absolute deviation from the mean normalised to the standard deviation is 0.79 and 0.90 for the PIC181 and PIC184 measurement respectively, and 2.76 and 1.73 for the PIC181 measurement superimposed with 1 nF and 1 mH. Defining a threshold value for this quantity can thus help to indicate artefacts in the measurement data that can not be explained by the model before the parameter identification process is started.

4 Material parameter identification

To efficiently solve the inverse problem of determining the piezoelectric material parameters from a given frequency-dependent, electrical impedance, gradient-based local optimisation is applied. A prerequisite for the application of local optimisation is the availability of sufficiently accurate initial values for the material parameters. The analytical approximations used to estimate initial values in previous work [8, 10, 11] are replaced by a neural network trained to approximate the inverse of the simulation model used to solve the inverse problem. Because the parameter estimated using the trained neural net constitute approximations, the classical optimisation procedure is still performed, however, with significantly reduced complexity and computational effort.

An overview schematic of the proposed optimisation strategy is given in Figure 5. A first estimate for the material parameters is determined using the inverted forward model (Sect. 4.1) with the measured impedance vector Z_meas as input. The resulting values p_init are then used as a starting point for two consecutive, gradient-based optimisation steps, that only differ in the forward model used: The first step employs the trained surrogate model (Sect. 4.2) and the second the slower but more accurate finite element (FE) model. Both steps iteratively modify the parameter vector p_i to gradually increase the similarity of the model output (Z_sur and Z_sim respectively) and the measurement data Z_meas.

Figure 5 Schematic of the developed two-stage optimisation strategy to determine piezoelectric material parameters p based on a measured electrical impedance vector Zmeas.

4.1 Parameter estimation network

The network for parameter estimation based on the electrical impedance is used to approximate the inverse of the simulation model (Sect. 2). The mathematical relation between input (the parameter vector p) and output (the electric impedance vector Z_sim) is highly non-linear and thus not classically invertible. Previous studies on initial value estimation methods using neural networks in the context of piezoelectric material characterisation indicate that a modified convolutional network architecture is suitable for this application [21]. The basic idea is that the convolutional layers, in which the input data is correlated with a number of trained filter kernels [24], are used to identify the shape of the local extrema resulting from resonance and anti-resonance effects. The typical convolutional network architecture includes pooling layer to reduce the networks cross-section. This pooling is typically performed in the dimension of the input data and thus leads to a translational invariance of the output of the network with respect to features in the input data. For the given application, however, the positions of resonance frequencies in the impedance data are prominent quantities influenced by the material parameters. Translational invariance should thus be avoided for the parameter estimation network. Pooling layers are, however, required to reduce the networks cross-section and thus implemented orthogonal to the dimension of the impedance data. The resulting network architecture (Fig. 6) employs two convolutional layers followed by maximum pooling. The logarithmic absolute value of the impedance as well as the phase are used as inputs for the network. The density ρ of the samples can be measured directly and is thus also used as an input, however, it bypasses the convolutional part of the network and is fed into the first of three fully connected layers. These layers gradually decrease the cross-section of the network. Only continuous, differentiable activation functions are used to promote smoothness of the inverted model. Hyperbolic tangent activation is used for the convolutional and dense layers. To not limit the range of the output of the model, the exponential linear unit is used for this layer. The outputs of the network are the material parameters scaled identically to previous approaches that employed classical optimisation strategies [7, 10].

Figure 6 Modified convolutional neural network architecture designed and trained to approximate the inverse of the simulation model.

Noise is added to the training data set for the training of the neural network to account for similar effects in the measurement data. The data set is split into 67% training and 33% validation data. Training is performed using adaptive momentum estimation [23] with a batch size of 256, a learning rate of 10⁻⁴, and a weight decay of 0.05 for 1000 epochs with the mean absolute percentage deviation as objective function. The training process (Fig. 7) shows successively decreasing costs for both training and validation data with no significant indication of overfitting occurring, which is observed when applying pure feed forward networks to a similar problem [21].

Figure 7 Value of the objective function (mean absolute percentage deviation) during the training process of the material parameter estimation network for the training and validation data.

The validation of the network’s performance for previously unknown data is possible because the exact inverse of the network exists in the form of the simulation model used to generate the training data. A new impedance data set is generated for arbitrary material parameters and supplied to the network. The resulting material parameters, which should ideally be identical to the previously assumed parameters, are then used to simulate the electrical impedance of the sample. A comparison between the originally simulated impedance and the impedance simulation with the material parameters estimated by the neural network shows good agreement (Fig. 8). Even smaller features and resonancefrequencies are reconstructed using the estimated material parameters.

Figure 8 Synthetic electrical impedance data generated using the simulation model for arbitrary material parameters compared to the impedance generated with the material parameters estimated by the neural network (Fig. 6).

Similarly, electrical impedance data from physical measurements Z_meas can be supplied to the network to estimate material parameters, given that it is sampled at the same frequencies. The resulting parameters p_init can be used to simulate a frequency-dependent, electrical impedance to assess the accuracy of the estimates (Fig. 9). The simulation using the estimates provided by the neural network shows adequate agreement with the measurement data, with the largest deviations observed around the thickness resonant mode at 2.1 MHz. Compared to initial value estimation based on evaluating analytical approximations for the frequency-dependent impedance [11], however, the neural network estimates show qualitatively similar results. A quantitative evaluation using the objective function applied for the subsequently implemented inverse procedure (Eq. (10)) results in significantly reduced deviations for the neural network estimates p_init (26.8) when compared to the analytic initial value estimation method ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ (117.5, cf. Tab. 2).

Figure 9 Measurement data for the electrical impedance of a piezoelectric ring (PIC184) compared to the impedance simulated with the material parameters estimated by the neural network pinit and using analytical estimates for the material parameters ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ [11]. Evaluating the objective function (Eq. (10)) shows that pinit yields an improved approximation of the physical data compared to ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ (values 26.8 vs. 117.5 respectively).

4.2 Surrogate model network

Aiming to further improve the performance of the optimisation procedure, a neural network is trained to mimic the behaviour of the finite element model of the piezoelectric sample. The goal here is to use the surrogate model for an intermediate optimisation stage before the finite element model is used in the final optimisation process. The surrogate model thus has to provide an adequate approximation of the finite element model. The computational effort to evaluate the surrogate model is, however, several orders of magnitude smaller than the evaluation of the finite element model. This intermediate step can thus be used to improve the estimates provided by the parameter estimation network without requiring significant computation time once the surrogate model network is trained. The training data can be used once more; the only large computational effort required for this intermediate step is thus the training of the model itself. The network used for this is a simple feed forward architecture using only fully connected layers. The input of the network are the material parameters including the density. Because the objective function for the subsequent optimisation procedure will be based on logarithmic absolute value of the impedance [8, 10, 21] the surrogate model is also trained to output only this, neglecting the phase. The model has five hidden layers with increasing numbers of neurons (200, 300, 500, 1200, and 1500), each again with a hyperbolic tangent activation function. The output layer has the size of the measured impedance vector (1999 neurons) with an exponential linear unit activation. The model is trained for 1000 epochs using adaptive momentum estimation [23] and the mean squared error as the objective function with a batch size of 256, a learning rate of 10⁻⁴, and 33% validation data.

Because the training of a neural net constitutes a gradient-based optimisation procedure, the respective machine learning libraries provide facilities to access the gradient of the network’s output with respect to the input values. This property of neural networks can be exploited in the subsequent gradient-based optimisation procedure, because the gradient of the surrogate model’s output Z_sur with respect to the input values does not have to be approximated, e.g. using finite differences, but can be extracted from the network. This further enhances the performance of this intermediate step, provided that the surrogate model network is a good approximation ofthe finite element simulation model.

A trust region reflective optimisation algorithm [25] is used to minimize the objective function (the squared difference of the logarithmic absolute value of the simulated and measured impedance vectors [8, 10, 21])

$$\begin{array}{c}\hfill \mathit{J}={({log}_{10}|{\mathit{Z}}_{\mathbf{meas}}|-{log}_{10}|{\mathit{Z}}_{\mathbf{sur}}|)}^2\end{array}$$(9)

starting from the estimates p_init provided by the parameter estimation network. The optimisation problem is bounded using the parameter ranges of the training data (Tab. 1), assuming the surrogate model is not able to extrapolate. Because the inference of the surrogate model is computationally inexpensive, this intermediate optimisation step can be performed in few seconds, yielding p_sur.

4.3 Optimisation with finite element model

The surrogate model only provides approximations of the physical behaviour of the sample based on training data and statistical processes. These models are as such to be considered unreliable to realistically describe the behaviour of the physical sample and only provide means to accelerate the quantitative material identification process. The finite element model, capable of accurately modelling the physical phenomena occurring in the sample, thus has to be evaluated at least once to verify the results provided by the machine learning-based methods presented here. It is, however, proposed to perform a final optimisation step using the finite element model of the sample. Previous work employed a block coordinate descent approach, partitioning the optimisation problem to subsets of the parameter vector before optimising the problem for the complete parametervector p [8, 10, 11]. A recent study, using a similar approach proposed here, shows, however, that the block coordinate descent method is not necessary if the initial values are of sufficient accuracy [21]. The optimisation is again performed applying a trust region reflective algorithm [25] to the objective function

$$\begin{array}{c}\hfill \mathit{J}={({log}_{10}|{\mathit{Z}}_{\mathbf{meas}}|-{log}_{10}|{\mathit{Z}}_{\mathbf{sim}}|)}^2\end{array}$$(10)

starting at the output values of the surrogate model optimisation p^sur and assuming no bounds. The output of this last optimisation step are the final parameter estimates p^res.

5 Results

The performance of the material parameter estimation method is evaluated by applying it to measurements of the electrical impedance of two samples made of different materials: PIC181 and PIC184 by PI Ceramic. Both materials are doped lead zirconate titanate ceramics designed for usage in high-power acoustic applications, so-called hard piezoelectric ceramics with low damping. Measurements are conducted using an impedance analyser (E4990A by Keysight) at frequencies matching the simulation data (1999 samples up to 8 MHz). As an example, the determined material parameters for the PIC181 sample at 25 °C at the different stages of the optimisation procedure (p_init, p_sur, and p_res) are listed in Table 2. The results ${\mathit{p}}_{\mathbf{res}}^{\mathbf{\prime }}$ when omitting the intermediate optimisation step and using the initial values ${\mathit{p}}_{\mathbf{init}}^{\mathbf{\prime }}$ for the gradient-based optimisation directly are also listed as an alternative method with reduced complexity. To compare the method with previous results, estimates based on analytical approximations ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ and the respective optimisation results ${\mathit{p}}_{\mathbf{res}}^{\mathbf{″}}$ are also shown [11]. The norm of the objective function (Eq. (10)) is evaluated using the finite element model to quantify the deviation between measurement data and simulations using each material parameter set.

A comparison between initial values determined from analytical approximations ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ and the estimates provided by the parameter estimation network p_init shows a significantly reduced objective function value, confirming the qualitative assessment based on Figure 9. Optimisation using the surrogate network as the model of the sample yields p_sur and further decreases the output value of the objective function. The procedure terminates after 27 iterations, each requiring a few milliseconds to execute on current standard desktop hardware. Of special note here is that the output for the mass-dependent damping parameter α_M is close to zero after these steps, even though the value of the objective function indicates that the physical behaviour of the sample is better represented by p_sur than by p_init. This behaviour confirms the results of the sensitivity analysis (Fig. 2): The sensitivity with respect to the damping parameters is very low. The second optimisation step using the finite element model, however, yields a non-zero result for α_M, with a further reduced objective function value. This result is produced after 20 iterations, each requiring 13 evaluations of the objective function to determine the gradient with respect to each material parameter along with the primal value. An evaluation of the finite element model(provided with the data set [17]) requires approximately one minute using current standard desktop hardware when parallelised in frequency domain. The complete parameter estimation procedure can thus be performed in few hours with no further parallelisation using the approach.

Alternatively, the initial values ${\mathit{p}}_{\mathbf{init}}^{\mathbf{\prime }}$, determined using the parameter estimation network, can directly be used as initial values for the gradient based optimisation procedure using the FE model, which now converges after 36 instead of 20 steps. The results of this simplified approach ${\mathit{p}}_{\mathbf{res}}^{\mathbf{\prime }}$ are remarkably similar to p_res, with significant deviations only occurring in the damping parameters α_M and α_K. That the objective function yields approximately the same value for p_res and ${\mathit{p}}_{\mathbf{res}}^{\mathbf{\prime }}$ again shows the low sensitivity of the procedure with respect to these parameters. This simplified strategy indicates, however, that the training and the optimisation with the surrogate model can be omitted at the cost of an increased number of evaluations of the FE model. It also indicates a certain robustness of the optimisation procedure with respect to the initial values.

A comparison with the results of the previously used method [11], that relies on analytical approximations for estimation of initials values ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ and block coordinate descent (BCD) steps in the gradient-based optimisation procedure, shows, however, that the final parameter estimates ${\mathit{p}}_{\mathbf{res}}^{\mathbf{″}}$ are highly dependent on the initial values. This indicates that the objective function has a number of local minima of different prominence. The approach also requires a significantly higher number of steps (321) in the optimisation procedure while yielding higher objective function values for the estimated parameters. This further highlights the requirement of adequate initial values for the piezoelectric material parameter identification problem.

To estimate the uncertainty of the identified material parameters, the finite element simulation model can be linearised at p_res using finite differences¹. Calculating a pseudo-inverse of the linearised model enables a means to propagate the uncertainty of the primary measurement quantities (the electrical impedance, the geometry and the mass of the sample) to the estimates for the material parameters [27]. The procedure can be adapted for the determination of piezoelectric parameter measurement uncertainty [18] and is compatible with the present study. The results of the uncertainty analysis for p_res are listed in Table 2. The relative standard uncertainties ${\mathit{u}}_{{\mathit{p}}_{\mathbf{res}}}$ are low for the elastic parameters $c_{11}^{\mathrm{E}}\dots c_{44}^{\mathrm{E}}$ as well as for the dielectric and piezoelectric parameters in 3-direction (${\epsilon }_{33}^{\mathrm{S}}$ and e₃₃). This observation aligns with the results of the sensitivity study (Fig. 2) and previous results for soft piezoelectric ceramics [8, 10]. Similarly, ${\epsilon }_{11}^{\mathrm{S}}$, e₁₅, and e₃₁ show increased uncertainties.Generally, the range of the standard uncertainties ${\mathit{u}}_{{\mathit{p}}_{\mathbf{res}}}$ of the non-damping parameters is higher compared to previous results [10] for piezoelectric material parameters. These previous results are, however, based on a sample geometry optimised for sensitivity, which is not the case for this study. The larger range of the values for the relative uncertainties is thus expected. Of special note are the uncertainties of the damping parameters α_M and α_K, which are very large. Especially the uncertainty for α_M indicates that this parameter can change in the order of a magnitude without significantly influencing the simulation and thus impedance result. α_M is, given the range typically observed for this parameter [8, 10], very low for the materials examined in this study. It can be assumed to be zero without influencing the objective function significantly. This would reduce the damping behaviour to the Kelvin-Voigt damping model [28], which is quantified by a single parameter, e.g. α_K. The uncertainty for α_K is a magnitude lower than for α_M, however, it is still considered large. In contrast to α_M, α_K can not be assumed equal to zero but has to be included to account for the present, but low, damping of the physical sample. The quasi-independence on α_M may also be explained be the fact that this quantity only influences the damping in the low frequency regime. As the sample is primarily modelled for elevated frequencies, the influence of α_M is further reduced.

That the uncertainty for α_K is lower than for α_M appears inconsistent with the results of the sensitivity analysis (Fig. 2), in which the impedance shows a larger sensitivity with respect to α_M than to α_K. This observation can, however, be explained by the larger, approximate value used for α_M to analyse the sensitivity. A larger value for the damping parameter results in larger changes to the overall behaviour of the virtual sample for given relative variations of this parameter.

A graphical comparison of the measured impedance and the output of the finite element model (Fig. 10) for p_res as well as ${\mathit{p}}_{\mathbf{res}}^{\mathbf{″}}$ shows good agreement with all major characteristics of the impedance accurately represented. The position of the local minima and maxima at the piezoceramic’s resonance and anti-resonance frequencies are identified accurately for both parameter sets. This shows that both parameter sets yield a suitable approximation of the physical behaviour of the sample, even though the numerical evaluation of the objective function indicates a better fit for p_res. The local extrema, however, are less pronounced in both simulation results compared to the measurement data, indicating that the model shows the same resonance frequencies but higher damping than the physical sample. This deviation may be attributed to the Rayleigh damping model’s inability to model the complete, broadband physical behaviour of the sample. The results already discussed show that the damping behaviour of the sample is primality identified by the α_K parameter, with α_M being negligibly small. The damping characteristics possible in the model are thus limited to isotropic effects that scale with frequency. Alternative damping models are discussed [10] and subject of future studies. Deviations observed at the resonance frequencies may also appear more significant than they are due to limited resolution in the depiction.

Figure 10 Comparison of impedances from measurements and from a simulation using the identified material parameters for a piezoceramic (PIC181) showing good agreement.

5.1 Temperature dependence

To assess the ability of the proposed method to detect small changes in material behaviour and perform reliable estimations for a larger number of measurements, it is applied to impedance measurements of the same sample at different temperatures. Impedances measurements are conducted in a temperature-controlled environment [26]. After allowing the samples to reach thermodynamic equilibrium, impedance measurement data is acquired from 25 °C to 85 °C in increments of 5 K. The material parameter identification process is applied to each measurement separately, yielding data sets for each sample and each temperature step. The gradient-based optimisation using the finite element model converges for each temperature step and both samples with objective function values of the same order of magnitude. The results for the PIC184 sample are visualised in Figure 11. Observations generally indicate a low influence of the temperature on the overall material behaviour. Especially the mechanical parameters $c_{11}^{\mathrm{E}}\dots c_{44}^{\mathrm{E}}$ show only little influence with respect to the changing temperature. Larger changes are observed in the permittivity, especially in ${\epsilon }_{33}^{\mathrm{S}}$. The results for ${\epsilon }_{11}^{\mathrm{S}}$ show erratic changes with respect to temperature, which can be attributed to statistical effects resulting from a low sensitivity (Fig. 2) or elevated uncertainty (Tab. 2) with respect to this parameter. Similar effects, although significantly less pronounced, occur in the piezoelectric constants e₃₁ and e₁₅ for the same reasons. The increase in e₃₃ with temperature can be considered significant, because of the distinct influence this parameter has on the sample’s behaviour. The damping parameters are both observed to decrease with temperature, but also show superimposed statistical deviations as expected given the results of the uncertainty analysis.

Figure 11 Results for the temperature dependence of the determined material parameters of PIC184 with linear regression for each parameter.

Even though the identified changes in the material parameter with respect to temperature are small, especially for the mechanical properties, they are relevant in the context of high-power acoustic applications. In these devices, the piezoelectric components are typically driven in a specific resonant mode, requiring the electronic hardware to be tuned to the resonance frequency. Even small changes, e.g. in the mechanical properties influence the resonance frequency, requiring frequency control hardware. To quantify these changes, the coefficients of linear polynomials

$$\begin{array}{c}\hfill p_i(\vartheta )=p_{i,25{}^{°}\mathrm{C}}+p_i^{\prime }·(\vartheta -25{}^{°}\mathrm{C})\end{array}$$(11)

are identified for each parameter p_i and both samples, where ϑ is the temperature. The results of this regression are listed in Table 3. A comparison between the two samples indicates that the mechanical quantities of PIC181 are marginally higher than those of PIC184, however, the temperature dependence is less pronounced in PIC181. The sign of the gradient $p_i^{\prime }$ is different for some quantities ($c_{13}^{\mathrm{E}}$ and $c_{44}^{\mathrm{E}}$), which may be interpreted as qualitatively different material behaviour between the samples. However, a more likely reason for this observation is that the effects caused by the temperature dependence are in the same magnitude as the uncertainty of the procedure with respect to this parameter and are thus stochastic in nature. An observation that aligns with the properties of the materials advertised by the manufacturer is the low damping observed for the PIC181 sample. While the damping generally decreases with temperature the parameter α_M is identified to be small and negative at 25 °C but increases with temperature and is positive starting at 27 °C. Comparing the primary results for α_M shows that the highest value observed for the PIC181 sample is 7.9 ms⁻¹ at high temperatures, which is significantly lower than the results for PIC184 (cf. Fig. 11). These results again indicate that the damping behaviour of the PIC181 sample can be accurately described assuming α_M to be zero.

6 Conclusions

The inverse measurement procedure supported by data-driven methods presented in this study allows for a robust and fast identification of piezoelectric material parameters. The resolution of the procedure is sufficient to show the temperature dependence of the material parameter in a moderate temperature range. Compared to works on similar problems [8, 10, 11], which require multi-stage optimisation with block coordinate descent procedures tailored to the specific sample configuration, the complexity of the identification procedure is significantly reduced. The computational effort required for the parameter identification, and thus the overall execution time of the optimisation procedure is also reduced compared to previous methods by an order of magnitude. However, because a training data set is required for the neural networks, the overall computation resources required by the procedure are significantly increased. In a production context, where parameters for the same device configuration are to be determined a number of times, this drawback may be negligible. The synthesis of the training data, however, is only feasible if sufficient computational resources, e.g. a computing cluster, are available. It should further be noted that by bypassing the analytical initial value estimation step using machine learning methods, which could rightfully be considered a brute-force approach, the knowledge typically acquired by the thorough physical and mathematical analysis of the problem may remain undiscovered.

The determined material parameter values and their temperature dependence can be used for physically accurate therm-piezoelectric simulation of the material’s behaviour. The quantitative results further indicate that for low-damping, hard piezoceramic materials like PIC181 the damping behaviour can be accurately described by the Kelvin-Voigt dampingmodel.

Future extensions of the method will include a variant of the identification procedure that is not limited to one sample geometry. This will require the synthesis of a larger training data set with varying geometry, however, similar to the density, the geometric parameters of the samples can be measured directly. Thus, the number of material parameters to be determined will not increase. Due to the different nature of the resonance behaviour, it is, however, still advisable to use different networks and optimisation strategies for different sample topologies, e.g. discs and rings. A study on reducing the amount of training data and finding the minimum number required for a successful model inversion within the context of piezoelectricity is also considered, as well as more advanced sampling methods. Future work will also include the identification of other damping models [10], which more accurately represent the physical damping behaviour over a wider frequency range. Further analysis will show if the approach presented here can be transferred to other inverse problems in measurement, for example acoustic waveguide-based approaches for polymeric samples [29] with time domain signals as measurement quantities. These studies will show if general guidelines for the machine learning-supported solution of inverse problem with application in measurement can be inferred.

Acknowledgments

The authors would like to thank the Paderborn Center for Parallel Computing (PC²) for the support and the access to the Noctua 1 and Noctua 2 clusters.

Funding

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 444955436 (research group 5208, NEPTUN).

Data availability statement

The training data sets [17, 19] and corresponding finite element models are available at https://zenodo.org/. The source code of the multi-physics simulation tool https://opencfs.org/ used in this study is available at https://gitlab.com/openCFS.

References

Appendix A Notes on uncertainty propagation

The uncertainty of the proposed method to determine piezoelectric material parameters can be estimated by linearising the simulation model and determining the pseudo-inverse [18, 27]. The parameter estimation network constitutes a similar concept to the pseudo-inverse of the simulation model, however, it is not limited to a linear mapping of the input data. The network is also not guaranteed show similar behaviour to the simulation model at arbitrarily small scales. However, because gradient information, e.g. the derivative of the output parameter p_i with respect to a sample Z_j of the impedance $\frac{\partial p_i}{\partial Z_j}$, is generally encoded and accessible for a trained neural network, the uncertainty of the input quantities (u(Z_j) and u(ρ)) can be propagated to the output of the network without significant computational effort. This approach can be used to estimate the influence the primary measurement’s uncertainty has on the final measurement quantities and may yield uncertainties for the initial value estimates. This approach will not account for the uncertainty of the network itself, which can be determined by evaluating training performance and is expected to be higher than the influence of the primary measurement quantities for the present study. Given the network’s shortcomings, however, this approach will yield propagated uncertainties with machine precision and does not rely on finite differences.

All Tables

All Figures

	Figure 1 Sectional view of the simulation model of a piezoelectric ring used in the present study. Due to rotational symmetry, a two-dimensional mesh of the cross-section area A can be used.
In the text

Figure 2 Results of a sensitivity analysis of the numerical model for piezoelectric discs with full electrodes (Full), the piezoelectric rings (Rings) studied here [17], and the piezoelectric discs with structured electrodes optimised for sensitivity (Structured) [8, 19]. Material parameters are varied by ±1% and the squared change in logarithmic absolute value of the electrical impedance is observed.

In the text

	Figure 3 Randomly selected examples for the synthetic impedance data from the data set [17].
In the text

	Figure 4 Distribution of the values in the latent variables of the trained autoencoder for the training data set (box plot) and values for specific measurements of different piezoelectric ceramics (PIC181 and PIC184, PI Ceramic) plain and with superimposed artificial measurement artefacts.
In the text

	Figure 5 Schematic of the developed two-stage optimisation strategy to determine piezoelectric material parameters p based on a measured electrical impedance vector Zmeas.
In the text

	Figure 6 Modified convolutional neural network architecture designed and trained to approximate the inverse of the simulation model.
In the text

	Figure 7 Value of the objective function (mean absolute percentage deviation) during the training process of the material parameter estimation network for the training and validation data.
In the text

	Figure 8 Synthetic electrical impedance data generated using the simulation model for arbitrary material parameters compared to the impedance generated with the material parameters estimated by the neural network (Fig. 6).
In the text

Figure 9 Measurement data for the electrical impedance of a piezoelectric ring (PIC184) compared to the impedance simulated with the material parameters estimated by the neural network pinit and using analytical estimates for the material parameters ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ [11]. Evaluating the objective function (Eq. (10)) shows that pinit yields an improved approximation of the physical data compared to ${\mathit{p}}_{\mathbf{ana}}^{\mathbf{″}}$ (values 26.8 vs. 117.5 respectively).

In the text

	Figure 10 Comparison of impedances from measurements and from a simulation using the identified material parameters for a piezoceramic (PIC181) showing good agreement.
In the text

	Figure 11 Results for the temperature dependence of the determined material parameters of PIC184 with linear regression for each parameter.
In the text