Parameter estimation in modelling frequency response of coupled systems using a stepwise approach

This paper studies the problem of parameter estimation in resonant, acoustic fluid-structure interaction problems over a wide frequency range. Problems with multiple resonances are known to be subjected to local minima, which represents a major challenge in the field of parameter identification. We propose a stepwise approach consisting in subdividing the frequency spectrum such that the solution to a low-frequency subproblem serves as the starting point for the immediately higher frequency range. In the current work, two different inversion frameworks are used. The first approach is a gradient-based deterministic procedure that seeks the model parameters by minimising a cost function in the least squares sense and the second approach is a Bayesian inversion framework. The latter provides a potential way to assess the validity of the least squares estimate. In addition, it presents several advantages by providing invaluable information on the uncertainty and correlation between the estimated parameters. The methodology is illustrated on synthetic measurements with known design variables and controlled noise levels. The model problem is deliberately kept simple to allow for extensive numerical experiments to be conducted in order to investigate the nature of the local minima in full spectrum analyses and to assess the potential of the proposed method to overcome these. Numerical experiments suggest that the proposed methods may present an efficient approach to find material parameters and their uncertainty estimates with acceptable accuracy.


Introduction
The scope of the current work is to investigate challenges arising in the solution of optimisation problems associated with the the performance of coupled dynamic systems or the estimation of their model parameters. Of particular interest is the occurrence of local minima, often posing a major obstacle in the solution of such problems.
There exist a vast amount of works published within the fields of material parameter estimation and optimisation. A comprehensive review of the literature is beyond the scope of the present paper; the interested reader is referred to [24,16]. Most works discussing performance optimisation in structural-acoustic systems either work on the eigenfrequencies or limit the cost function to frequency bands containing a few resonances, see e.g. [28,14,29,6,1,30,13,21].
However, the problem of interest here, i.e. estimation of model parameters in resonant, coupled problems over a frequency range containing a number of resonances, remains open. In particular those cases where, over the range of frequencies studied, there is a varying degree of interaction between the physical fields involved, are widely agreed to be subjected to local minima.
These often renders the parameter identification questionable if at all tractable, as demonstrated in previous works by the authors and others, see e.g. [23,4,15,5,20,12]. The presence of resonances and anti-resonances in the spectrum is not the only factor responsible for local minima. In fact, it has been observed in previous works, see for example [15,5], that a sufficient degree of anisotropy is prone to induce such phenomena as well. The underlying cause for this is thought to be the self-similarity of the structural response at discrete combinations of the material, geometrical and loading properties. However, in order to distinguish between local minimas related to such dynamic response phenomena, and those related to the resonant spectrum fitting in itself, robust and efficient combinations of inversion methods are interesting to study.
The present paper addresses a method for overcoming local minima in the estimation of material properties in coupled dynamic resonant problems over a wide frequency band. The proposed approach relies on a stepwise solution, where the solution of each sub-problem is the starting point of the next. While a multi-step concept as such is not new, see for example [27] where a 3-step procedure was employed to enhance convergence, [10] where different models are employed to fit sub-sets of parameters, and [19] where seismic waves are used to estimate subsurface properties using a frequency-by-frequency inversion approach, the proposed strategy is novel in that it specifically addresses resonant problems by employing an incremental frequency domain augmentation approach as discussed in the paper.
In our work, we have chosen to combine two different estimation methods, one based on a gradient, local search approach, [22], and one global search based on the Bayesian inversion framework, [24,11,2]. This enables us to combine the best of two analysis frameworks, [3,20], complementing the deterministic search with additional insight into ill-posedness, measurement, and/or model uncertainties, etc. in the inverse problem. This adds a major advantage as the probability density computed from the Bayesian solution, reveals not only the different point estimates but also tells about the correlation between different parameters and also the credible intervals of the estimates. In addition, while the Monte Carlo based estimate converges to the global minimum if the number of samples is not limited, the combination provides more accurate initial conditions for the chain(s) and this holds the potential of reducing the computational time.
The proposed method is here demonstrated for a simplified one-dimensional fluid-structure problem exhibiting a complex behaviour, chosen for its computational tractability as well as being inexpensive to solve. The studied problem possesses a known solution and its fluidstructure coupling effects are limited to the interface between the two media. This provides the possibility to conduct extensive numerical experiments with controlled characteristics, such as frequency resolution, dissipative losses, nature of coupling, etc. The model problem chosen for the current work is selected based on the authors' previous experiences from estimation in resonant systems involving soft materials, and the difficulties associated with local minima [5,25,26]. In the latter works, parameter estimation has required manual guidance to a certain extent in order for the cost function to converge. The authors aim herein at assessing the convergence behaviour of the proposed methodology, thus opening its applicability to more complex cases.
The paper is structured as follows. The model problem formulation is given next in Section 2. The deterministic and Bayesian inversion frameworks are introduced in Section 3 while Section 4 is dedicated to numerical experiments. Finally, discussion and conclusions are given in Sections 5 and 6, respectively.

Background and problem formulation 2.1 Model problem
Consider a freely moving one-dimensional (1D) solid body, excited at one end by a harmonic force of amplitude F , and coupled to an air cavity which is rigidly terminated at its other end, see Fig. 1. Both domains are modelled in the form of the 1D wave equation and have the same cross-section area A. Let p denote the acoustic pressure in the air cavity, u the linear solid displacement, ρ s the solid density, and ρ f the fluid density. The domain lengths are set to L s and L f for the solid and the fluid, respectively. The harmonic time dependence is chosen as e +jωt , where ω is the angular frequency. The solid and the fluid are both assumed to be damped, with loss factors η s and η f , respectively. The complex speed of sound for the solid and the fluid may thus be written as whereĉ f,s denotes the amplitude of the speeds of sound. The complex wave numbers for the fluid and solid are then The solution to the coupled problem is straightforward and will not be detailed here. The mean sound intensity level, at frequency ω, radiated from the structure to the fluid is given by where I ref = 10 −12 W·m −2 . The sound intensity level I is used in the current work as the target quantity that the model should reproduce as a basis for the parameter estimation. In Eq. (3), the pressure and the solid displacement (at x = 0) are given by In Eqs. (3)(4)(5), ξ denotes the set of known parameters and θ the set of unknown parameters,

Target properties
As stated in the introduction, the problem chosen is selected on the basis of it being simple to model and having a low computational cost, yet complex enough to exhibit local minima, which pose a formidable challenge to most approaches for inverse modelling. The fixed parameters for the solid have been chosen to be representative of a fairly light open-cell polymer foam, with a density low enough to induce a strong coupling with the air cavity over a wide frequency range.
The parameters for the target model are shown in Table 1, where θ meas denotes the target parameters governing the measured system. Table 1: Fixed parameters ξ and target parameters θ for the studied problem.
Fluid Solid Global With the chosen target model parameters, the fundamental frequencies and first harmonics of the uncoupled resonances of the fluid and the solid are shown in Table 2. Note that the fundamental frequencies of the two systems are well separated, and hence implying a moderately strong interaction between the fluid and the structure for these. Also, the fundamental resonance of the fluid cavity is within 1 Hz from the first harmonic of the solid, which leads to the wellknown up-and down-ward shifts in the coupled frequencies occurring for every second resonance of the solid, as will be shown below. As the focus is here on the phenomena rather than computation efficiency, a frequency range of 0.5 − 1000 Hz is used, with a frequency resolution of 0.0153 Hz, in the numerical experiments performed.
3 Inverse problem

Observation model
Recalling the aim of this work, together with the choice of model problem as described above, we consider also the influence of noise in the target spectrum data, from now on referred to as measurement data. This choice is partially made in order to avoid the "inverse crime" trap discussed below, but also to render the measurement data set more realistic and in this way obtain a preliminary assessment of the robustness of the proposed method. For the purpose of the current work, the noise is introduced in the form of additive errors on the sound intensity. The observation model may thus be written as where ω (ν) ∈ R N f ×1 (N f denotes the number of frequencies) and e models the measurement error. In this work, error e is assumed to be Gaussian distributed with zero mean and covariance Γ e = σ 2 e I, that is e ∼ N (0, Γ e ). In the covariance matrix, σ e denotes the standard deviation of the noise and I is the N f × N f identity matrix.
As stated in the introduction, we have chosen to apply and combine two different approaches to the inverse problem, a gradient search method and a method based on a Bayesian framework. Both of these are well-known and have been discussed in numerous papers, thus we will not go into details about either of them as the focus here is on the estimation and the problem of local minima in resonant systems.
Each method is run several times in order to get some statistics of their respective convergences to the target model. For each such repeated computation, the starting point is generated from a uniform distribution between the lower and the upper bounds of each of the variables respectively.

Numerical target
In the current work, a purely simulation-based generation of the measurement data is used and these are fitted with an analytical model. In such studies, it is of great importance to avoid committing an "inverse crime" [11]. Our choice for avoiding this is to use a discretized solution, leading to a different numerical accuracy for generating simulated measured data and the model fitted in the inversion. In this paper, we thus generate the measured data set using the Comsol Multiphysics software version 5.3a while in the inversion we use the analytic solution of the problem.
For the Comsol model, which consist of an axi-symmetric model, the element size is L e r = L e x = 0.0058 m for the solid elements in a mapped quadrilateral mesh with second order shape functions, while for the fluid we choose elements with L e x = 0.0344 m and L e r = 0.0058 m. Thus the mesh is compatible along the interface between the fluid and the solid. For both media these element lengths observe L e x,r < λ min /10 for λ min = c s /f max . With increasing frequency, the relative quadratic error between the discretized Comsol and the analytical solutions grows to around 7%, as shown in Fig. 2.

Deterministic framework
For the gradient search, the Globally Convergent Method of Moving Asymptotes (GCMMA) is used in the current work [22]. The cost function to be minimised is defined as where θ ∈ A, with A = [θ min , θ max ]. The parameters θ are not subjected to any additional constraints. We employ a re-scaling of all design variables such that they are mapped to the dimensionless domain θ ∈ [1, 2], see [4] for more details. The gradients required for the cost function are computed using finite differencing, with a step of dθ = 10 −12 .
In the analysis, a termination criterion of the GCMMA iterations is set to 10 −8 for the variation of the cost function and 10 −6 for the change in the design variables, as monitored over 5 subsequent iterations. As the focus is on the problem of local minima, and not on computational efficiency, the number of iterations for each gradient solution is allowed to reach 600 in order to ensure that the iterations are terminated based on convergence, to facilitate the evaluations described below.
In the following, the solutions obtained from the gradient search method will be referred to as GCMMA.

Bayesian framework
In the Bayesian framework [24,11,2], estimated variables θ are modelled as random variables. In this inference, the main tasks are to construct the likelihood and prior models π(I meas | θ) and π(θ), respectively. The likelihood model gives the relative probabilities that the model generates the observations I meas over all choices of parameters, while the prior model sets constraints on the model parameters. The solution to the inverse problem in the Bayesian framework is coded in the posterior distribution π(θ | I meas ), which gives the relative probability of the unknowns θ given the observations I meas . The posterior distribution is obtained by Bayes' formula π (θ | I meas ) ∝ π (I meas | θ) π(θ).
The additive observation model (8) allows us to write the likelihood model as follows where π e is the probability density of the error e. Since the error e is assumed to be Gaussian, the likelihood density can be written as Furthermore, we use uniform prior of the form: For the sampling of the posterior densities, we use the Metropolis-Hastings algorithm [18,17,9] with an adaptive proposal distribution scheme [7,8]. The adaptive algorithm automatically adjusts the spatial orientation and size of the proposal density considering all of the target distribution points accumulated so far, leading to a faster convergence.
As an initial condition for the Markov chain, we use either the GCMMA solution or random starting point and run it for 1,030,000 iterations including a burn-in of 30,000 (unless stated otherwise in the text). From the burn-in, the last 20,000 samples are used to provide an initial estimate of the covariance matrix for the adaptive proposal distribution scheme. As the proposal density, in the first 30,000 samples, we use a zero mean Gaussian distribution, in which parameters are assumed to be uncorrelated.
The point estimate for the parameters we use in this work is the maximum a posteriori (MAP) which corresponds to the point in the posterior density that has the highest probability. As an indication of the uncertainty in the point estimates we use the (narrowest) 95% credible interval, which represents the spread of the posterior density. In the following, the solutions obtained from the Bayesian inversion will be referred to as MCMC (Markov chain Monte Carlo).

Numerical experiments
This section presents the results obtained for the local minima problems in focus here. First the inverse problem for the whole frequency range is solved, this will be referred to as the full spectrum estimation, including an analysis of the character of the local minima as such. The proposed stepwise method is then introduced as a means to overcome these problems.
In the following examples, estimated parameters are restricted to the interval, η s ∈ ]0, 1], η f ∈ ]0, 1], c f ∈ [10, 700] m/s, and c s ∈ [10, 700] m/s which is referred to as the prior model in the Bayesian framework. This choice is made in order to avoid adding a priori knowledge with respect to each of the coupled sub-systems while still keeping a reasonable range of possible parameter values.

Full spectrum estimations
In order to explore the nature of possible local minima resulting from GCMMA, we attempted an estimation of the parameters given in Table 1, for a noise-free target spectrum generated through a Comsol simulation as described above. Figure 3 shows the initial and final parameter values for 392 different full spectrum estimations using random starting points, none of which converged to the global minimum. The speed of sound for the fluid tends to be close to the target value for most of the solutions, while the speed of sound for the solid seems to converge towards values that are considerably higher than the target value. These results illustrate the well-known difficulties encountered when fitting a parametric model over a wide range of frequencies, containing a large number of coupled resonances. In fact, in order to identify the parameters in Table 1 on the full spectrum, the initial values for θ must lie in the close neighbourhood of the target values.

Investigating the local minima
In this section, the focus is on two randomly chosen samples from Fig. 3, with the intention to study the results from the full spectrum estimations in more detail. We have selected starting points that are fairly close to the target value for the fluid speed of sound, but far away from the solid, the first having close to 10 times higher and the second nearly half the target value. Thus, the points The results are shown in Fig. 4 in terms of sound intensity level. For both θ 1 and θ 2 ,ĉ f is close to the value used in the target model, while for neither of the two design pointsĉ s is close to the target model value. In addition, as observed in the close-up y ∈ [0, 0.01] in Fig. 3, none of the estimates for the solid loss factor have converged to the target value, see Table 1.
The solutions obtained are quite typical for the two different inversion approaches. With the restrictions set, primarily in terms of number of iterations, GCMMA terminates at different local minima which clearly are far from the true solution, while MCMC identifies the target model but often needs a long chain to find the global minimum. Note that in Fig. 4, the MCMC solver does not use adaptivity. For the purpose of this investigation of the local minima, the proposal density is assumed to be wide, allowing the chain to jump between different local minima more easily.
In order to confirm whether the obtained GCMMA solutions correspond to local minima of Eq. (9), the vicinity of θ full 1 and θ full 2 is explored using MCMC. Here, the MCMC algorithm uses the adaptive proposal density after the 30,000 initial rounds but it is assumed that the proposal density in these initial rounds is narrow in order to restrict θ to the neighbourhood of these GCMMA solutions. Figures 5 and 6 show that both solutions obtained are indeed local minima. All four design variables influence the cost function to the same extent and obviously there is no direction in which a gradient search can expect a lowering of the cost function further. In addition, it can be observed from these two figures that the correlation between the model parameters vary between almost totally uncorrelated (see e.g. in Fig. 5ĉ f andĉ s ) and higher correlation (e.g.ĉ f and η f ; η f and η s ). Furthermore, the 95% credible intervals reveal that the local minima are indeed narrow in all parameters.
It is worthwhile to explore the solutions obtained for the full spectrum in some more detail. In almost all cases we have observed, the first resonance in the response spectrum (see Fig. 4) is well identified with respect to the frequency but not particularly well in terms of amplitude. The gradient solutions appear in many cases to track the uncoupled resonances, see for example the peak in the GCMMA solutions in Fig. 4 which occurs at about 172.4 Hz and the peak at 342.2 Hz, both very close to the uncoupled fluid resonance frequencies given in Table 2. In most of the solutions, this appears to be achieved through converging towards design points where, e.g.ĉ f ≈ĉ target f /2 orĉ f ≈ĉ target f ; andĉ s ≈ 2, 4, 5, 8 ×ĉ target s , see Fig. 3. An interpretation of these solutions is that they tend to combine resonances in one of the uncoupled systems with anti-resonances in the other. Together with associated high loss factors this leads to the fitting to the anti-resonant parts of the target spectrum that may be observed in the lower part of Fig. 4.
As a concluding remark of this subsection, similar findings have been observed for a case with target loss factors of 10%, and these will not be repeated here for the sake of conciseness.

Stepwise approach
A major point of attractiveness of GCMMA for solving inverse problems is its low number of model evaluations [22]. However, as illustrated above, a pattern of local minima hinders the unsupervised determination of the global minimum or the full spectrum problem. Observations made when studying full spectrum solutions such as those in Sec. 4.1 triggered the idea for the stepwise approach presented here for solving inverse problems involving frequency-dependent experimental data exhibiting a resonant behaviour. The method draws on the observations made related to the correlation between the model parameters as well as the convergence of the gradient method with respect to the lowest resonances of the coupled system. The concept is simple, yet, as detailed in what follows, remarkably powerful and robust.
The first step of the method consists in defining a number of sub-problems in which the problem is solved within a frequency interval from the lowest available frequency to an incrementally In the very first step, i.e. when solving for F sub 1 , θ init is randomised such thatĉ s andĉ f are sampled with equal probability from the interval [10, 700] m/s and damping factors η s and η f from the interval ]0, 1], i.e. the prior model. In each of the subsequent -th steps, θ init holds the parameters estimated from the solution obtained in the previous step, i.e. θ init = θ −1 . That is, the point estimate obtained from the first subdivision is used as a starting point for the next frequency interval until the entire available spectrum is covered. The last step corresponds to an estimation in the full spectrum F sub N sub where the starting point has been guided by the sequence of previous steps.

Application to the simplified fluid-structure interaction problem
This subsection presents the application of the proposed stepwise approach to estimating the parameters of the target model given in Table 1. The focus is on the convergence of the method and therefore the computational efficiency aspects are not discussed at this point.
The application of the proposed approach to the dataset shown in Fig. 4 yields systematically satisfactory results, only limited by the numerical stopping criteria set for the optimiser. As these results do not convey a significant deal of information, the performance and robustness of the method are presented on a dataset containing additive noise. This is done in order to demonstrate how the stepwise approach performs for a more realistic set of data, by including a standard deviation σ e = 2 dB SIL (where SIL stands for sound intensity level) assumed for the zero mean white noise model e in Eq. (8).
In setting up the estimation problem, the frequency spectrum f ∈ [0.5, 1000] Hz is divided into N sub = 5 frequency intervals, having maximum frequencies f max = {5, 25, 75, 150, 1000} Hz, see Eq. (14). This choice, which obviously is far from unique, is based on the observations made in Sec. 4.2 on the lowest peak in the response occurring at circa 17.5 Hz and the correlations betweenĉ f andĉ s . In particular, the first substep is aimed at capturing the low frequency asymptotic behaviour and the parameters governing this, in this caseĉ f and η f as will be discussed later.
Using θ init = θ init 2 , given in Sec. 4.2, as a starting point, the results in Fig. 7 are obtained. The fit between the target and the estimated model is very close, and for comparison, the less satisfactory estimation results when using the full spectrum for the inversion are also shown. Figure 7: I as a function of frequency: model including additive noise, final step in the proposed approach, and full spectrum estimates following Table 3.
The convergence history for the stepwise solution is shown in Table 3, together with the full spectrum estimates and initial and target parameters. The stepwise solution and the target values are also visualized in Fig. 8. Interestingly the path to the target solution is by first finding the loss factor for the fluid while the loss factor for the solid goes to the maximum allowed. By the third step, all remaining parameters have been estimated. Note that this is achieved for a quite noisy target dataset. In analysing the convergence behaviour of the stepwise approach, it should be kept in mind that in each substep the subproblems solved are different from the full spectrum problem. Nevertheless, the solutions to the incremental steps provide a path of incremental partial global minima towards the global minimum of the full spectrum problem.  Table 3. Horizontal black lines show the target parameter values.  Fig. 9, indicate that there is a strong correlation between the fluid loss factor and the fluid speed of sound, but the remaining parameters are neither correlated to each other nor influencing the solution to the subproblem as such.

Marginal densities for the different steps
At the 3th step, the correlation is strong betweenĉ s andĉ f , and it can be observed that both have indeed converged to the target values. As the full spectrum is solved in the final step, all parameters have been estimated within 5% relative error, with the largest deviations occurring in the loss factors, in this case related to the noise in the measured data set used for the estimation. The figures also show how the credible interval becomes narrower when the amount of data is increased for all parameters. The true value is found within the narrowest 95% credible intervals in all other cases, except in Fig. 11 forĉ s but also here the distance between the point estimate and the true value is acceptable, given the uncertainties in the measurement data.
To study the robustness of the stepwise approach when applied to the model problem studied in the current work, we performed 392 runs, using the same randomly chosen starting points as in Fig. 3. Out of these 9 did not converge to the global minimum, see Fig. 12. In the converged solutions, the speeds of sound for both fluid and solid, as well as the loss factor for the fluid are estimated within 1% relative prediction error, while in particular the loss factor for the solid is less accurately estimated. The largest discrepancy observed for the latter is about 15%.

Discussion
As mentioned above, the success of the proposed stepwise approach appears to be linked to the dominant model parameters and the correlations between these and their variations over the different parts of the frequency spectrum studied. This is in agreement with the findings in [27,10].
In fact the correlation between parameters observed above for the current model problem, may be illustrated further by taking the limit for the Helmholtz numbers k f L f , k s L s 1, of Eq. (8). This may be shown to be Accordingly, as observed above, for the first substep in the solution of the current model problem, c f and η f are the two model parameters that control the behaviour of I in Eq. (3). This in itself motivates our choice of the first sub-divisions in the previous section, and in particular Figure 10: Results for the the substep f ≤ 75 Hz, otherwise same caption as in Fig. 9.
the choice of the frequency range for the first substep, which in many of the converged solutions is critical for the success. By restricting the first substep solution to essentially fitting Eq. (15) to the target spectrum, the ratio between the fluid loss factor and the square of the fluid speed of sound, which according to Eq. (15) is governing the intensity at low frequencies, is obtained, i.e. ∼ 4.237610 −09 , already at the end of the first step as ∼ 4.157310 −09 . As a matter of fact such ratio is invariant along the path in the marginal density forĉ f vs. η f , Fig. 9. A crucial aspect for the interpretation of the solution to the first step is the fact that the marginal densities related to the solid loss factor are nearly uniform. This indicates that the cost function is not sensitive to this parameter for the first step, which in turn suggests that the sound intensity radiated by the solid into the air cavity does not strongly depend on η s at very low frequencies, as confirmed by Eq. (15).
In addition, Eq. (15) and Fig. 9 suggest that the low-frequency asymptotic behaviour of the system becomes increasingly sensitive to decreasing values of the speeds of sound. This translates into the difficulty for a gradient-based optimiser to approach the correct solution from substantially low values of the speeds of sound. Figure 11: Results for the last substep f ≤ 1000 Hz, otherwise same caption as in Fig. 9.
As identified in the convergence of repeated runs from random starting conditions, Fig. 12, the estimates related to the solid loss factors are the least accurate. We relate this to the numerical model in Comsol, used as measurement data, and the computational errors associated with it, see Fig. 2 where it is shown that the solid displacement is the dominating source of error for most frequencies in the spectrum. Together with the added noise, this renders the convergence in the final substep mostly to be controlled by the solid loss factor, see for example Table 3 and Fig. 11.
The choice of the lower bounds of the parameter range used, in particular for the speeds of sound, was dictated by our observations when studying the full spectrum convergence behaviours. In fact, the robustness of the approach relies on the ability of the first step to fit the asymptotic behaviour. Lower speeds of sound would accordingly require a lower high-frequency limit for the first step.
Although its general applicability is still to be confirmed by exploring real parameter estimation cases, the significant improvement in convergence experienced for the current coupled fluid-structure problem is indisputable in our opinion. We ran almost 400 analyses with uni- Figure 12: Left: Sample starting points that did not converge. Right: Relative prediction error for the converged samples. The solid parameters are denoted by "×" and the fluid by "•".
formly distributed random initial starting points, all but a few converged to the target solution.
Even though the focus of the current work is not on computational efficiency, it is worth noting that as the frequency range of the first steps is relatively narrow, it is far less expensive to compute the cost function for these as compared to the full spectrum approach. Thus, it should be beneficial to allow for a higher number of iterations, i.e. employing more stringent convergence criteria in the first increments of the stepwise estimation, in order to reduce the efforts required in the last steps. This is obviously a trade-off that most likely is problemdependent as well as depending on the starting point and the path required to reach the global minimum.

Conclusions
Solving estimation problems over a wide frequency range with several resonances is difficult due to local minima, which arise due to partial fits to the spectrum. Such an obstacle is particularly severe for gradient search methods. Nonetheless, also global search approaches face difficulties, requiring a large number of iterations.
In the present paper, the proposed stepwise inversion approach is run using several hundreds of random starting points in the parameter design space. The method is observed to converge to the target solution in the vast majority of cases. This represents a substantial improvement of the gradient-based parameter search on its robustness to a systematic pattern of local minima due to a resonant behaviour. In addition, although the convergence rate depends on the initial starting point and hence the computational cost required to obtain the solution, the stepwise approach proposed in this paper provides potentially a means to reduce the number of model evaluations required to achieve convergence, by reaching increasingly good convergence as the evaluation frequency is increased.
Our work in the current paper used a combination of MCMC and GCMMA in order to understand the convergence behaviour. The settings applied for each of the two methods were partly dictated by this. Nevertheless, in the authors' opinion, the central findings in this paper are independent from the optimisation tool at hand. The stepwise method proposed herein provides a means towards circumventing local minima in full spectrum inverse estimation problems involving resonant structures by exploiting their asymptotic behaviour at low frequencies.