Introduction

Currently, real-time measurement of parameters within an electromagnetic field, such as internal heat generation or heat capacity, is technically unfeasible. To overcome this situation, those parameters can be estimated by solving the corresponding inverse problem. Besides, it is now possible to work out these type of problems using diverse global optimization algorithms (^{Hào et al., 2017}; ^{Wang et al., 2017}; ^{Cebo-Rudnicka et al., 2016}; ^{Nedin et al., 2016}; ^{Chen et al., 2016}; ^{Strongin & Sergeyev, 2000}; ^{Zhigljavsky & Žilinskas, 2008}; ^{Kvasov & Sergeyev, 2014}). Previous works related to heat generation proposed estimating heat generation in the case of a rotatory friction welding system (^{Yang et al., 2011}). To estimate the time-dependent heat generation at the interface of cylindrical bars during the rotary friction welding process, they used an inverse algorithm based on the conjugate gradient method and the discrepancy principle. ^{Bermeo et al. (2015)} focused on a different application, i.e. heat cancer treatment, and they demonstrated that by using radiofrequency electromagnetic waves, it is possible to solve the corresponding inverse problem and estimate state variables, such as temperature distribution in the tissues. The authors arrived at the solution using the Sampling Importance Resampling (SIR) algorithm. Recently, and dealing with a heat conduction process, Wang *et al.* (Wang & Liu, 2016) estimated thermal conductivity using a nonlinear parabolic equation with a temperature-dependent source. They solved the inverse problem using an iterative optimization algorithm. Similarly, ^{Tutcuoglu et al. (et al.2016)} estimated the thermal parameters of internal Joule heaters. Their results were experimentally verified by using a block of polydimethylsiloxane embedded with a strip of conductive propylene-based elastomer. The inverse scheme suggested is based on the governing nonlinear, inhomogeneous heat conduction and generation equation. ^{Mohebbi et al. (2016)}, using the conventional *conjugate gradient method* and the two-dimensional inverse heat conduction problem, estimated thermal conductivity, heat transfer coefficient, and heat flux in irregular bodies. Equally, they described a sensitivity analysis for the solution of the inverse problem. An inverse problem, oriented to the identification of the heat capacity and thermal conductivity, was solved by Nedin *et al.* (2016). They solved it using an iterative algorithm for the Fredholm's equations of the first kind. ^{Hussein et al. (2014)} in their paper explain the inverse problem of determining the time-dependent thermal conductivity from Cauchy data, considering one-dimensional heat equation with space-dependent heat capacity. The model, a parabolic partial differential equation, was solved using the standard finite-difference method. On the other hand, Bermeo *et al.* (2015) solved an inverse problem using several PSO strategies. Their goal was to estimate the particle size distribution of colloids from multiangle dynamic light scattering measurements. Additionally, ^{Giraldo et al. (2012)} used the inverse problem in order to estimate neural activity from electro-encephalographic signals.

More recently, ^{Mohebbi et al. (2016)} proposed how to estimate parameters such as thermal conductivity, heat transfer coefficient, and heat flux in three-dimensional irregular bodies in steady state for a heat transfer conduction process. They used a sensitivity analysis scheme for computing the corresponding sensitivity coefficients in a gradient-based optimization method. Similarly, they asserted some advantages when using this sensitivity analysis, such as simplicity and accuracy, which makes the solution of the inverse problem very accurate and efficient. Their solution strategy started by generating an elliptic grid, then solving Fourier steady state heat equation using the traditional finite-difference method to determine the temperature value at each node. This article also included a diversity of references published in the last decades and dealing with the assessment of such thermal parameters. Likewise, ^{Cui et al. (2016)} modified the traditional Levenberg-Marquardt method using a complex-variable-differentiation approach to do the sensitivity analysis. This new version seems to have the same advantages of the original algorithm, but it also yields better convergence stability and efficiency due to its accurate evaluation of the sensitivity coefficients. As observed from the literature review, it is possible to establish that inverse heat transfer problems, dealing with the estimation of thermodynamic materials properties, is still a common and useful strategy.

Sometimes, as in the present case, it results unfeasible to measure properties such as heat capacity, due to factors that affect these experimental determinations. Measuring temperature and magnetic or electric field intensities, for example, in an electromagnetic environment, is quite tricky today. In the present article, we propose a strategy for the estimation of parameters such as heat capacity and volumetric heat generation during heat treatment of a material. The manuscript includes a brief description of the algorithms required for the solution of the inverse problem, along with the description and solution of the direct problem. After that, some of the more relevant results are presented and analyzed. At the end, we include the most relevant conclusions.

Materials and methods

This section includes a brief description of the traditional Levenberg-Marquart (LM) optimization algorithm, and modern optimization algorithms such as Spiral Optimization Algorithm (SOA), Vortex Search (VS) and Weighted Attraction (WAM). We also propose hybrid algorithms (SOA-LM, VS-LM, and WAM-LM) that synthesize the best advantages of the traditional and modern approaches. Moreover, the direct and inverse problems are stated.

*Algorithm fundamentals*

The *LM* method is a well-known iterative deterministic technique traditionally used for non-linear parameter estimation. Details are omitted for the sake of brevity, but they can be found elsewhere, e.g. in Necati ^{Ozisik & Orlande (2000)}. In contrast, *SOA* (^{Tamura & Yasuda, 2011}), *VS* (^{Dogan & Olmez, 2015}), and *WAM* (^{Friedl & Kuczmann, 2015}) are considered global optimization metaheuristic algorithms that try to mimic natural phenomena such as pressure fronts, vortex pattern created by the vertical flow of the stirred fluids, and gravitational attraction between particles in order to solve optimization problems. It is important to emphasize that one can distinguish between the global optimization from the local one, because the first one focuses on finding the extreme of a function in the whole search space (function's domain) (^{Kvasov & Mukhametzhanov, 2017}; ^{Sergeyev & Kvasov, 2017}). Below, we briefly describe these algorithms.

Spiral Optimization Algorithm (SOA)

SOA consists of five steps:

Step 0: Algorithm initialization. Establish the number of spirals in the solution space, the number of maximum iterations, the rotation angle and the convergence radius.

Step 1: Spirals center selection. Evaluate the initial point of each spiral in the objective function (OF). Then, choose the one with minimum value as the rotation center.

Step 2: Spirals rotation. Rotate the remaining spirals around the spiral center selected in the previous step.

Step 3: New spirals center selection. Evaluate the new set of points in the OF, and choose the one with the minimum value as the new center.

Step 4: Stop criteria check. If the convergence criteria are satisfied, the algorithm stops. Otherwise, it goes back to Step 2.

Vortex Search (VS)

VS consists of four steps:

Step 0: Algorithm initialization. Establish the number of particles in the solution space, the number of maximum iterations, the initial center and the initial radius. The last two parameters are based on the limits of your search space.

Step 1: Distribution of the particles. Generate possible solutions by using Gaussian distribution around the center with a standard deviation (radius).

Step 2: New center and radius selection. Evaluate possible solutions in the OF. Then select the best solution to replace the current center. And then, decrease the radius for the next iteration according to the gamma distribution.

Step 3: Stop criteria check. If the convergence criteria are satisfied, the algorithm stops. Otherwise, it goes back to Step 1.

Weighted Attraction Method (WAM)

WAM consists of five steps:

Step 0: Algorithm initialization. Establish the number of particles in the search space and the number of maximum iterations and explosions. Then, randomly place every particle in the search space. Besides, set as zero initial movement distance.

Step 1: Evaluate OF. Evaluate candidate solutions in the OF and assign an attraction factor to each one of them.

Step 2: Movement of the particles. Calculate the center of the mass of the particles and then move the particles to that center based on their previous and current movement.

Step 3: Explosion. If particles are too close to each other make an explosion (disperse the particles randomly again).

Step 4: Stop criteria check. If the convergence criteria are satisfied, the algorithm stops. Otherwise, it goes back to Step 1.

On the other hand, we propose new hybrid methods between the aforementioned classical algorithm and the metaheuristics. The main idea of these hybrid methods is to first run a metaheuristic, and then use its solution as the initial conditions of the LM method. The basic flowchart of this method is shown in Error! Reference source not found. With these hybrids, we expect to preserve the stochastic behavior of solutions from metaheuristics while maintaining the proved convergence of the classical method.

*Mathematical model of the system*

For demonstrative purposes, an isotropic and homogeneous solid sphere of radius (a), with constant density *(p)* and specific heat (c) is considered. The heat conduction equation expressed in spherical polar coordinates assumes no variation for the (θ, *φ* ) coordinates. The rate of internal heat generation per unit volume at *r = 0* is *q(r")* . The mathematical model is, thus, given by 1, where *k* is the thermal diffusivity of the substance; *K* is the thermal conductivity of the substance, and *T (r,t)is* the temperature. The relationship between *k* and *K* is given by *k* = *K/p.* Temperature at its boundary and in the initial condition is assumed to be 25°C.

Statement of the direct problem

It is possible to calculate the temperature distribution within a sphere homogeneously irradiated by microwaves. Here, we assume that all parameters of the above mathematical model are known, including the initial and boundary conditions. In the same way, we select three functional forms for the heat generation parameter, i.e. linear, sinusoidal and exponential.

Statement of the inverse problem

Based on the solution of the direct problem, it is possible to estimate properties such as heat capacity as well as internal volumetric heat generation parameter *q'*
_{
Q
}
*"* within the solid sphere. Solving this problem requires experimental temperature readings at various radial positions within the sphere, and knowledge of the nature of the internal heat generation . Nonetheless, they can be gathered using one or several temperature sensors.

Direct problem solution

This problem already has a well-known general analytical solution given in (^{Carslaw & Jaeger 1959}) that yields the temperature profile within the sphere as shown in Equation (2), where for the linear form, for the sinusoidal form, and for the exponential form. The term q_{0} is the magnitude of the internal volumetric heat generation term. In order to plot this solution, Equation (2) was normalized for the variables *(T, r* and, *t)* and the new ones are denoted by *(T*
_{
u
}
*, r*
_{
u
} and, *t*
_{
u
}
*)* as shown in Equation (3), where each term is detailed in Equation (4).

Figure 2 shows the normalized theoretical temperature within the sphere as a function of the normalized radius and normalized time, as well as the profile temperature after a long time of irradiation. As expected, the temperature field is axisymmetric, having the maximum value at the center. For the figure, numbers on the curves represent the parameter *t*
_{
u
}
*.* For the sake of brevity, the other two graphs for the sinusoidal and exponential case were omitted.

Simulating temperature measurements

We used synthetic temperature values for simulating real measured values, with the objective of checking the validity of this approach for a microwave heating process analysis. As expected, real temperature measurements may contain random errors. Here, such errors are assumed to be additive, uncorrelated, and normally distributed with zero mean and a constant standard deviation. Hence, we accept as valid the basic standard statistical assumptions proposed by Beck (^{Vere Beck & Kenneth 1977}). Consequently, a white Gaussian noise (AWGN) was added to the theoretical normalized temperature to construct synthetic ones. In this case, SNR was kept constant at 30 dB. Figure 3 shows an example of the theoretical and measured temperatures for a given position and time. In the figure, temperature measurements at normalized time tu = 0,1, are shown as an example. For conciseness, the other two graphs for the sinusoidal and exponential case were omitted.

Objective function (OF)

The OF that provides the minimum variance in an inverse problem is the ordinary least squares (OLS) norm. In this work, we are going to take into consideration the case when the transient readings (Y) were taken from multiple sensors, as shown in Equation (5), where Y and T are vectors containing the synthetic (measured) and theoretical temperatures, respectively. It is worth remarking that due to the simulated nature of this work, whenever temperatures are labeled as "measured" they actually indicate simulated measurements generated by our model.

Error and RMS error

We considered throughout this work two metrics for the error. The reasoning behind this decision was that sometimes algorithms may converge to an estimated solution distant to the theoretical one , but for which the evaluation of the objective function (OF) is similar. Moreover, by using both metrics we can explain why in some cases the theoretical and estimated parameters exhibit a high percentage of error in the parameters, while having a low RMS error in the evaluation of the OF. Therefore, the first one was the error of the parameters themselves, as shown in equation 6, where theoretical parameters are contrasted against the estimated ones (i.e. ״ and cp). The second one was the RMS error of the temperature profile, as shown in equation 7, where the theoretical temperature (T) is contrasted against the estimated temperature (Test). The latter was calculated with the parameters obtained through the optimization algorithms. In other words, the RMS error relates, in a way, the evaluation of the OF1 to OF2.

Results and analysis

This section is divided into two main subsections. The first one describes the application discussed above. Within this, the results obtained with the algorithms executed on the data (i.e. LM method, metaheuristic algorithms, and the hybrid algorithms) are presented. The second one is a performance analysis that was carried out with the results derived from all the algorithms.

A demonstrative example

Consider the microwave treatment of an isotropic and homogeneous silicon carbide (SiC) sphere with constant density, thermal conductivity, and specific heat respectively (Grup d'Innovació per la Millora de la Docência en Estructura Propietats i Processat de Materials n.d.). Its diameter is 2 [cm] and internal heat generation rate at *r* = 0 is 8000 [W/m^{3}]. The temperature at its boundary is fixed at 25 °C. Furthermore, it is assumed that four sensors are located at unitary radius *r*
_{
u
} = 0,3, *r*
_{
u
} = 0,5, *r*
_{
u
} = 0,7 and *r*
_{
u
} = 1. The goal is to estimate *c* and , by solving the inverse problem. It is important to keep in mind that most of the reported data are normalized. Therefore, parameters c and obtained through the algorithm are also normalized. However, data shown within this section relate to their unnormalized equivalents.

LM method

Table 1 summarizes the results obtained with this algorithm. However, it is important to note that this algorithm converges under certain initial conditions (IC). If these IC are very far from the true values (i.e. 0.9 times the true value for the lower limit and 35 times for the higher limit in the linear case), LM does not converge. Table 2 shows the ranges tested for this particular problem.

Metaheuristic algotithms

This time, SOA, VS and WAM were used to find parameters *c*
_{
p
} and . To do so, several tests were carried out. The range of the IC where the metaheuristic algorithms were tested is shown in Table 3. The parameters used in each algorithm are shown in Table 4. Each algorithm was repeated 50 times and resulting data are summarized in Table 5.

**Hybrid algorithms**

In this scenario, the hybrid algorithms were used, considering the same parameters in Table 3 and Table 4. Tests were repeated 50 times and results were summarized in Table 6. As data show, using heuristics to identify proper initial points for *LM* leads to improved results (Table 2 and Table 5). Thus, these hybrids are more robust than the single use of heuristics.

**Algorithms performance**

A performance comparison of all algorithms is carried out. Results are shown in Table 7, Table 8 and Table 9. In these tables, it can be seen that LM algorithm requires less time and iterations than the metaheuristic and hybrid algorithms. Besides, hybrid algorithms require, in average, 10 seconds and 19 iterations more than metaheuristics. Additionally, Table 10 and Table 11 show an analysis of the accuracy based on the results of all experiments (i.e. 900). It is important to keep in mind that a guess is considered right when the parameter estimation resides within a 10% confidence interval. If it is outside, the guess is deemed wrong. Through these tables, we show that the results obtained with the hybrid algorithms improve the results obtained with the metaheuristics in over 60%.

Conclusions

This article described the estimation of the internal volumetric heat generation and the heat calorific capacity parameters during microwave heating of a solid spherical sample. We used three functional forms (linear, sinusoidal, and exponential) for the electromagnetic field. We incorporated White Gaussian noise into the direct problem to simulate a more realistic situation. Then, we solved the inverse problem through three different kinds of approaches. The first one used a traditional deterministic Levenberg-Marquardt *(LM)* algorithm. The second one used modern stochastic metaheuristics *(SOA-VS-WAM).* The final one was the proposed hybrid of the first two. Results show that *LM* converges to the expected solutions only if the initial conditions (IC) comply with specific ranges. However, LM only required 0,6% of the time employed by metaheuristics and hybrids. Nonetheless, through metaheuristics we can expand the range of the IC by 42 times. However, these methods only yield an accuracy of about 36%. Thus, by using hybrid algorithms we can merge the best of both approaches, increasing the accuracy of metaheuristics by about 60%, while only requiring about 6% more time. Therefore, we highly recommend to use this kind of hybrids for these inverse problems.