Genetic Algorithm as the Solution of Non-Linear Inverse Heat Conduction Problems: A Novel Sequential Approach

Direct measurement of surface heat flux in engineering applications is often a challenging if not infeasible, task. Examples include surface heat flux measurement for space vehicles during atmospheric entry, machining tools during operation, and thermal management of electronic circuits. The temperature measurement, on the other hand, can be performed effectively at subsurface locations. Heat conduction problems can be categorized as either direct problems or inverse problems. A direct problem involves using the known active boundary condition (e.g., heat flux, temperature) and solving the heat diffusion equation to find the temperature distribution with respect to time and space. The direct problem is mathematically well-posed and various solution methods for solving the resultant partial differential equation have been established. Such methods include the Green's function, separation of variables and Laplace transform, to name a few [1,2]. In the case of the inverse heat conduction problem (IHCP), the problem involves determining the unknown boundary condition using the known temperature from the interior points of a heat-conducting body. Such a problem is mathematically ill-posed and extremely sensitive to errors in measurement [3].

IHCP solution techniques can be considered in two separate categories: the whole-time domain method and the sequential method. The whole-time domain method considers the entire time domain to generate the solution and thus, this method can only be performed after an experiment has been undertaken. Examples of specific IHCP solution techniques that use the whole-time domain method include the Tikhonov regularization method, conjugate gradient method, and singular value decomposition [3]. The whole-time domain method is a possible solution approach; however, it is not useful for monitoring or control purposes where computing surface heat flux continuously in a near real-time fashion is needed. The other solution method, the sequential method, fills this role. An example solution method using the sequential frame of reference, known as the function specification method, calculated the surface heat flux using temperature data from past and future time steps, allowing the capability for continuous, near-real-time application. Following this route, it was shown that using the sequential method with Tikhonov regularization, the surface heat flux solution can be written in digital filter form [4], whereby the surface heat flux can be calculated by multiplying a filter matrix into a vector of measured temperature values. The authors Woodbury and Beck investigated the characteristics of the resultant filter matrix which showed that the rows of the filter matrix are all identical, only shifted in time, and filter values at both ends of the filter row approaches zero. The primary strength of the filter matrix produced by the digital filter method is its high computational efficiency; paired with it being a sequential method, it allows for the solving of IHCPs in a near real-time fashion requiring only a “window” of temperature data consisting of some previous and future time steps. Following its introduction, the digital filter method has been further developed by Najafi et al. for multilayer and multidimensional problems, as well as domains with temperature-dependent material properties [3,5–8]. Uyanna et al. developed the algorithm for an IHCP with a moving boundary [9,10], while Zhang et al. utilized an artificial neural network (ANN) to solve the IHCP with a moving boundary [11] and ablation [12].

When the problem is framed as such, the entire domain is fed into an IA which attempts to predict the boundary condition. The predicted boundary condition is then fed into the direct problem and a resultant temperature is computed and compared with the measured temperature. Subsequent iterations attempt to produce a boundary condition value that minimizes the error between the computed and measured temperatures until convergence. Once that is done, an optimal boundary condition time series is provided thereby solving the IHCP, such as the approach undertaken by Raudensky et al. [13] to find the unknown heat transfer coefficient subject to a boundary and Pourgholi et al. to find the unknown heat flux [14] both of which used the genetic algorithm (GA) for optimization. Likewise, using the ant colony optimization (ACO) algorithm and artificial bee colony (ABC) algorithm, boundary condition reconstruction was performed through similar problem constructions [15,16]. A review of some of these algorithms was investigated and presented by Udayraj et al. [17] who compared the performance of ACO, particle swarm optimization (PSO), and the cuckoo search algorithm, as well as Hetmaniok et al. [18] who compared ACO and ABC algorithms. Another popular solution method includes the machine learning algorithm known as the artificial neural network (ANN or NN), which can be applied using the whole-time domain approach as shown by Krejsa et al. [19], Raudensky et al. [20], and Szénási et al. [21]. A similar approach using a modified version of the ANN is presented by Deng et al. [22].

Genetic algorithms, inspired by natural evolution, have received a lot of attention for performing optimization in engineering applications, given their outstanding capabilities for finding global optima for nonlinear, multivariable, and constrained problems [23]. Particularly, genetic algorithms have been used for numerous heat transfer problems including heat exchanger optimization [24,25], designing thermal management systems [26], parameter estimation [27], and solving inverse heat transfer problems. The inverse identification of temperature-dependent thermal conductivity using a genetic algorithm was performed by Czél and Gróf [28]. A modified genetic algorithm for solving the inverse heat transfer problem for estimating heat source is studied by Liu [29]. Liu [30] also explored the application of particle swarm optimization for solving IHCP in order to calculate unknown surface heat flux in a circular channel. Pourgholi et al. [14] developed a numerical approach by combining the least squares method and the genetic algorithm for estimating the boundary temperature in an IHCP. A comprehensive review of the application of genetic algorithms for solving heat transfer problems is presented by Gosselin et al. [31].

The nature of the filter coefficient approach, which basically allows the calculation of surface heat flux using a limited time window of temperature values and a set of fixed filter coefficients, provides the opportunity of implementing it through machine learning algorithms. Najafi and Woodbury [32] used ANN in the form of a filter coefficient method and used it for solving IHCP in a 1D domain in a sequential manner. They used a known set of temperature (inputs) and heat flux (output) data and trained an MLP feedforward neural network, where the weights and biases basically could serve as the equivalent of filter coefficients. The trained network was then successfully tested for predicting the surface heat flux for several test cases using temperature data from a limited time window as the input.

In the present paper, the use of genetic algorithms as a method for the calculation of filter coefficients for solving IHCP in a one-dimensional domain is explored. A one-dimensional slab model was created in COMSOL Multiphysics for numerical simulation. Three heat flux profiles are applied to the surface of the model and temperature distributions are calculated in the domain (direct problem). The results from the solution of the direct problem are then used to calculate the filter coefficients using a genetic algorithm through minimization of the sum of the square of the errors between known heat flux values and calculated heat flux values using the filter coefficients. Three additional test cases are then generated and used to assess the performance of the proposed method in terms of accuracy, computational time, and robustness. The results show that the genetic algorithm can be successfully implemented for calculating filter coefficients to be used for solving one-dimensional IHCPs for domains with both constant material properties and temperature-dependent material properties.

Consider the schematic of the direct problem and inverse problem shown in Figs. 1(a) and 1(b), respectively. Transient heat flux is applied to the surface of the one-dimensional domain (x = 0) and the backside of the domain is insulated (x = L).

where 0 < x < L and t > 0. The initial and boundary conditions for this model are listed as follows:

The direct problem is solved using COMSOL Multiphysics by applying several heating profiles to the surface and calculating the temperature distribution within the 1D domain. The temperature data are then used to calculate the surface heat flux through the proposed inverse algorithm which will be discussed as follows.

where m_p is the number of significant filter coefficients before the current time-step and m_f refers to the number of significant filter coefficients after the current time-step. In other words, m_p and m_f are the number of temperature data from previous and future time steps in order to determine the value of surface heat flux at the current time-step (⁠$q̂M$⁠). The characteristics of the filter matrix are better illustrated in Fig. 2, which shows several of the rows of a filter matrix [Eq. (10)]. These filter coefficients are calculated using Tikhonov regularization method [Eq. (9)] using midplane temperature (x = L/2), a dimensionless time-step of 0.05, and Tikhonov regularization parameter of 0.01. As seen, the rows of the filter matrix (i.e., each curve on Fig. 2) are identical but shifted in time. Additionally, there are only limited number of significant (nonzero) filter coefficients and in each row. The components of the middle row of the matrix are further specified by markers in Fig. 2. Note that there are only about 20 nonzero filter coefficients, suggesting that calculating the heat flux at a given time-step, requires temperature data from 10 previous time steps (m_p=10) and 10 future time steps (m_f =10).

The primary strength of the filter matrix produced by the digital filter method is its high computational efficiency paired with its sequential nature, allowing solving IHCPs in a continuous near real-time fashion and requiring only a window of temperature data consisting of some previous and future time steps. It also allows addressing nonlinear IHCPs with temperature-dependent material properties through interpolating a series of filter coefficients as discussed by Beck [33], Najafi et al. [8], and Keltner et al. [34].

The use filter coefficient method requires calculating the sensitivity matrix, X [Eq. (9)], and proper selection of the Tikhonov regularization parameter. As discussed by Woodbury and Beck [4], one needs to know the nature/type of the heating profile (smooth changes, sharp changes, etc.) in order to properly select the order of magnitude of the Tikhonov regularization parameter, leading to accurate performance of heat flux calculation using filter coefficients. This might introduce challenges when there is no knowledge of the type of heat flux profile. The use of the filter coefficient method for nonlinear problems, although very powerful and computationally efficient, requires the knowledge of temperature-dependent material properties. This study explores whether a set of measured temperature-heat flux data can be used to find the optimal filter matrix through genetic algorithm to be used for solving 1D-IHCP and whether the solution will be robust enough to be applied to various types of heat flux profiles. This will be considered a significant enhancement to the filter coefficient technique as it can be implemented solely using temperature-heat flux data from one experiment, without the need for a full analytical solution to the direct problem, and applied to other cases as will be discussed later.

where $mp$ and $mf$ are the number of past and future time steps considered, respectively, $fJ$ is the filter coefficient pertaining to the temperature window index $j$⁠, and $Y$ is the temperature value at the time-step $M$⁠.

The filter coefficients are considered as optimization variables here. The filter coefficients are first calculated by using a known set of input (temperature) and output (heat flux) data and then the obtained filter coefficients are used to evaluate an unknown surface heat flux just by introducing the input data and using the previously calculated filter coefficients. This allows us to explore whether the filter coefficient calculated for a particular heat flux profile, will be capable of calculating surface heat flux when exposed to a different heat flux profile. This approach is comparable with NN training and testing procedures, where the network is first trained by a set of known input–output data and the weights and biases are determined accordingly [35]. The trained network can then be used as a tool for calculating outputs associated with a new set of inputs (predicting tool) [32].

Genetic algorithm is selected as the optimization method given its excellent characteristics for solving multivariable minimization problems with nonlinear objective functions as well as the superior performance in finding global optima [36]. Genetic algorithms are evolutionary algorithms inspired by natural evolution and has been successfully used for many different engineering problems [23,37], including heat transfer problems [31]. The optimization process starts with an initial population, which consists of a series of chromosomes (i.e., candidate solutions), and each chromosome consists of a group of genes (i.e., m_p+m_f filter coefficient values). The objective function is evaluated for each chromosome from the initial population and those with superior performance are identified and selected. In this case, the objective function is given in Eq. (12). A new generation is then created through reproduction, which involves crossover as well as mutation. During crossover, genes or subsets of the best members are recombined to produce new members of the population. The children of the best members are made of the genes of the best members, which will ideally replace the worst members in the population. During mutation, random changes occur to the genes in individual chromosomes to possibly improve the results. The objective function is then evaluated for the new generation. The process continues until a termination criterion is achieved.

It should be noted that in using a genetic algorithm for finding filter coefficients, it was not clear whether the filter matrix and filter coefficients were going to demonstrate similar characteristics as those calculated using the Tikhonov regularization method (e.g., see Fig. 2). In other words, the answers to three research questions need to be explored:

(1) Do the filter matrix and filter coefficients generated through the genetic algorithm demonstrate a similar pattern to those produced through Tikhonov regularization method? (2) Is the genetic algorithm able to calculate filter coefficients that can allow surface heat flux estimation in a near real-time fashion? Note that here “near real-time” refers to the use of temperature measurement data from internal surfaces from a limited time window (limited number of previous and future time steps) resulting in an acceptable delay. (3) Can filter coefficients that are evaluated from one given set of data (including temperature and heat flux) be used along with temperature measurement data to calculate surface heat flux associated with a different heat flux profile (i.e., assessment of the robustness)? The answers to these questions will be explored in the Secs. 3.3, 4.1 and 4.2.

It should be noted that for the optimization problem, the filter coefficients are calculated for each time-step to minimize the error between predicted heat flux and known heat flux at any given time-step. In other words, assuming values for m_p and m_f, a total of m_p+m_f coefficients are considered as the optimization variables to minimize the objective function given in Eq. (12). These coefficients are calculated for every time-step resulting in a filter matrix with N rows (N is the number of total time steps) and m_p+m_f columns.

A one-dimensional model was created in COMSOL Multiphysics in order to develop test cases for the purpose of this study. The physical properties of the model are given in Table 1. It should be noted that the material selected for this test study is Carbon Phenolic.

The temperature-dependent material properties used for the model are presented in Fig. 3.

A series of known transient heat flux profiles are then applied to the surface of the model and temperature values at the midplane are calculated and recorded. Three heat flux profiles (cases 1–3) were used for training. In other words, the heat flux values and corresponding midplane temperatures were used to find the optimal filter coefficients, which will result in the lowest value of the objective function, as defined in Eq. (12). Note that the number of filter coefficients (i.e., the number of optimization variables) is equivalent to m_p+m_f that covers the time window for which the current surface heat flux is impacted by the temperature values. This value is obtained through multiple trials with the goal of finding the lowest possible number of filter coefficients that lead to accurate results. In other words, the optimization is performed multiple times to evaluate filter coefficients assuming different values for m_p and m_f and the lowest number that resulted in good result were found. Notice that the smaller value for m_f is desirable as it will translate into a shorter delay in evaluating the surface heat flux. The heat flux and corresponding temperature values for cases 1–3 are shown in Figs. 4–6, respectively.

Three additional test cases are also developed using more complex shapes for the heat flux profiles. These heat flux profiles and their corresponding midpoint temperature values are shown in Figs. 7–9. The optimal filter coefficients are evaluated using the data for cases 1–3. Each of these sets of filter coefficients is then used to calculate surface heat fluxes for cases 1–6 (with added noise) using temperature values for each case and Eq. (13). Therefore, cases 1–3 are considered for training purposes (i.e., calculating optimal filter coefficients) and cases 4–6 are for cross-validation to explore whether the optimal filter coefficients can accurately estimate surface heat flux for different types of profile shapes. The heat flux profiles and the corresponding midplane temperature values for cases 4–6 are demonstrated in Figs. 7–9, respectively.

Upon calculating the filter coefficients using the genetic algorithm, the coefficients are plotted for each time-step (i.e., rows of the filter matrix), to explore the answer to the first question. Figures 10 and 11 show the selected rows of the filter matrix obtained for case 1 and 2, respectively. Here, the 50th, 150th, 250th, 350th, 450th, and 550th rows of the filter matrix are demonstrated. As seen, unlike the filter coefficients calculated using Tikhonov regularization, although the rows of the filter matrix look somewhat similar, they are substantially different. In addition, the filter coefficients do not approach zero from both ends. This indicates that the components of one row (one set of filter coefficients) cannot be used to calculate the surface heat flux for the entire time domain, unlike the case with the filter form Tikhonov regularization technique. Therefore, here the filter coefficients for each time-step are specific to that time-step. In other words, instead of using one set of filter coefficients for calculating the surface heat flux for the entire time domain, a filter matrix will be used with m_p+m_f columns and N-m_f rows, where the ith row of the matrix will be used. Here, N is the total number of time steps. Although this will translate into more computational cost since a larger number of filter coefficients must be calculated, since this process is done prior to the use of the solution, it should not impact it's usability.

The input heat flux profiles were used on the 1D slab modeled in COMSOL Multiphysics. As seen, test cases 1 and 4 both involve smooth variations while test cases 2 and 6 both involve sharp changes, and test cases 3 and 6 involve sudden jumps and constant heating. The physical properties used for the simulation are listed in Table 1.

It should be noted that for training and testing purposes, the calculated temperature values are modified with an addition of 0.25% random noise (representing measurement error) in postprocessing. This should present some difficulty to the algorithms' predictions and assess if the proposed solution method is robust against temperature measurement errors. It should be noted that the amount of noise is set based on the actual expected measurement error that may be expected when using a thermocouple (around ±2 K).

The results of all test cases and detailed discussions for each case are presented in this section under two subsections.

The data obtained from cases 1–3 were used to calculate the filter coefficients through the genetic algorithm. Three sets of filter coefficients are obtained accordingly (associated with cases 1–3). Temperature data (with added noise) were then used for cases 1–6 and filter coefficients were applied in order to calculate the surface heat flux. The parameters used for the genetic algorithm were found through trial and error to yield the best result. In other words, several combinations of parameters are used to find the best performance (highest R² value). These parameters are listed in Table 2.

The filter coefficients calculated using case 1 are used to evaluate the heat flux profiles associated with cases 1 and 4, to explore the answer to the second research question. This is because cases 1 and 4 are variants of a similar type of heat flux profiles (smooth variation). In other words, filter coefficients that are evaluated using the data from case 1 (training) are expected to perform well in predicating the surface heat flux for case 4 (testing). Similarly, the filter coefficients that are found through case 2 (training) are used to evaluate the surface heat flux for cases 2 and 5 (testing). Finally, filter coefficients that are found using data from case 3 are used to predict the surface heat flux for cases 3 and 6 (testing). The predicted heat fluxes for cases 1–6 are shown in Fig. 12. As seen, the heat flux estimation is quite accurate for all cases.

The coefficient of determination (R² value) between the predicted heat fluxes and actual values of heat flux are calculated and summarized in Table 3 to provide insights into the performance of the developed method. As seen high R² values (close to 1) are obtained for all cases which indicates that the predicted heat flux values are very close to the actual heat fluxes. Also, it can be seen that the highest R² value is determined for case 2 and the lowest one was found for case 3. It should be noted that the coefficient of determination is an indication of how close the predicted heat flux values are to the actual values.

Another study was performed to investigate the ability of the proposed approach to predict profiles which are not variants of the one that it was trained on. This is to respond the third research question. In other words, for the Robustness Test case, the optimal filter coefficients that are obtained on the simple parabolic profile (case 1) are then used to predict the other profiles (Cases 2, 3, 5, and 6). This was done to investigate if the algorithm has strength in generalizing the data that it was trained on for application toward other unknown profiles. Such an ability would be rather useful in real-world applications as there can be many different factors, environmental or otherwise, which can produce a heating profile that was not considered when training the algorithm. The results of this robustness test for Cases 2, 3, 5, and, 6 are presented in Fig. 13.

As seen in Fig. 13, even though a dissimilar heat flux profile (case 1-parabolic) is used for optimizing the filter coefficients, the predicted heat flux values for all other four cases (2, 3, 5, and 6) are very accurate which is an indication of the robustness of the proposed method.

The R² values associated with the prediction versus actual heat flux values related to the robustness tests were found as 0.9908, 0.9318 and 0.9788, and 0.9533, for cases 2, 3, 5, and 6, respectively.

It appears that the lowest R² values both for the testing and robustness test were associated with case 3 (train step variant). The highest R² value was found to be associated with case 2 (train triangular variant). The step function is known to be a more challenging profile to predict given the sudden jump that it experiences in the value of heat flux. As seen in Figs. 13(c) and 13(d), largest deviations are observed where a sudden jump is experienced in the step function. Note that the filter coefficients are evaluated using a parabolic heat flux profile and therefore may not result in very accurate results with a step function profile. However, the overall, the prediction looks very promising. The presented approach shows merit over Tikhonov regularization-based filter coefficient method in terms of robustness as the optimized filter coefficients can be obtained through one set of experimental data without the need for knowing the form of the heat flux profile. The interpolation of filter coefficients to account for temperature-dependent material properties is also no longer necessary as this information will be captured in the optimal filter coefficients. However, a matrix of filter coefficients must be optimized ahead of time in order to effectively use this method. Unlike the conventional digital filter method, the filter coefficients found by the genetic algorithm are not similar for different time steps. In terms of accuracy, it is expected that if used for a similar variant of a given heat flux profile, the Tikhonov regularization method will perform better; however, GA will provide more accurate results for other type of heat flux profiles and therefore considered more robust.

In the present paper, a novel approach is developed for using genetic algorithms to determine filter coefficients as the solution of IHCP in one-dimensional domains with temperature-dependent material properties. Six test cases were developed in COMSOL Multiphysics and the direct problems were solved through numerical simulation to calculate temperature distribution within the domain. Temperature data from the midpoint were used as the input and the surface heat flux was used as the output and the filter coefficients were optimized in order to minimize the error between actual heat flux values and the predicted heat fluxes. The optimized filter coefficients were then tested with other test cases both with similar and dissimilar natures of heat flux profiles to assess the performance and robustness of the proposed technique. The important findings are summarized as below:

• The characteristics of the optimal filter coefficients calculated by the genetic algorithm were compared against the coefficients that can be determined through the digital filter form of the Tikhonov regularization method. It was found that, unlike the conventional digital filter coefficients, the rows of the optimal filter matrix are not identical.

• The number of filter coefficients required for calculating surface heat flux at each time step was found as 50 (m_p = m_f = 25). This translates into a 25 s delay in calculating the surface heat flux. If smaller time delay is desirable, one possibility is to explore whether the sensor can be moved closer to the surface. This will allow using a smaller value of m_f and m_p.

• The filter coefficients that were calculated through the use of genetic algorithm applied on a set of heat flux-temperature data for a parabolic heat flux profile were successfully used for estimating various types of surface heat flux profiles in a near real-time fashion and with good accuracy.

An effective use of the developed solution is to calculate the filter coefficients for a given system with a set of known heat flux-temperature data. The filter coefficients can be then used to evaluate other unknown surface heat fluxes using internal temperature measurement for the same system. The presented approach offers a computationally efficient, accurate, and robust method for solving nonlinear IHCPs in one-dimensional domains in a near real-time fashion and can be applied for surface heat flux estimation and monitoring in a variety of applications including thermal protection systems of space vehicles, machining tools in manufacturing processes, and control of wall heat flux in furnaces and combustion chambers.