Estimation of space-dependent single scattering albedo in a radiative transfer problem using differential evolution

Abstract

This work presents the differential evolution algorithm applied to an inverse radiative transfer problem formulated as a finite dimensional optimization problem. It is considered a one-dimensional isotropically-scattering medium with finite optical thickness, space-dependent scattering albedo and plane-parallel geometry. The direct radiative transfer problem models the transmission of radiation through this medium by a linear version of the Boltzmann equation with azimuthal symmetry. The direct radiative transfer problem solution, which is required for the optimization techniques, is obtained by using the collocation method. Some test cases are presented, aiming at illustrating the efficiency of the methodologies used in the treatment of an inverse problem of radiative transfer. The results are compared with other approaches and indicate that the proposed methodology characterizes a promising methodology for dealing with this type of inverse problem.

Keywords:

• space-dependent scattering albedo
• radiative transfer
• inverse problem
• differential evolution

1. Introduction

The formulation and solution of direct and inverse radiative transfer problems are directly related to several relevant applications in a large number of areas of scientific and technological interest such as tomography 1,2, remote sensing and environmental sciences 3,4, and radiative properties estimation 5–9, radiative property identification of semi-transparent scattering media 10, identification of thermal radiation properties of dispersed media 9, to mention only a few.

Traditionally, inverse problems have been treated by using either classical or deterministic optimization techniques 11,12. In recent years however, the use of non-deterministic techniques or the coupling of these techniques with classical approaches thus forming a hybrid methodology became very popular due to the simplicity and robustness of evolutionary techniques 6,7,13–15. The main goal of this article is to solve the inverse radiative transfer problem with space-dependent single scattering albedo by using the differential evolution (DE) algorithm. To solve the direct problem the collocation method is used. The results obtained are compared with those obtained with other approaches proposed in the literature. This work presents a brief review on the DE technique. Both formulations for the direct radiative transfer problem and the associated inverse problem are presented. Afterwards, the methodology used in this work is described and applied to four illustrative examples. Finally, the conclusion and suggestions for future work are outlined.

2. Mathematical formulation and solution of the radiative transfer problem

In this work, a one-dimensional, grey, heterogeneous, isotropically scattering participating medium of optical thickness τ_o and transparent boundaries surfaces, as shown in Figure 1, is considered 16,17.

Figure 1. Schematical representation of a one-dimensional participating medium subjected to the incidence of radiation originated at external sources.

The boundaries at τ=0 and τ=τ_o are subjected to the incidence of radiation originated at external sources with intensities A₁ and A₂, respectively. Y represents the intensity of the radiation that comes out from the medium and may be measured by external detectors. The mathematical model for the interaction of the radiation with the participating medium is given by the linear version of the Boltzmann equation 16,18, which for the case of azimuthal symmetry and a space dependent albedo is written in the dimensionless form as (1) (2) (3) where I(τ, μ) represents the radiation intensity, τ is the optical variable, μ is the cosine of the polar angle, i.e. the angle formed between the radiation beam and the positive τ axis. The other symbols have already been defined, except for the single scattering space-dependent albedo, ω(τ).

When the geometry, the boundary conditions and the radiative properties are known, the problem given by Equations (1–3) may be solved and the radiation intensity I(τ, μ) is determined for the whole spatial and angular domains, i.e. 0 ≤ τ ≤ τ_o and −1 ≤ μ ≤ 1. This is the so-called direct problem.

To solve the direct problem, the collocation method 19,20 was used. In this methodology, the general Boundary Value Problem is described as (4) where x is the independent variable, y is a vector of dependent variables and p is a vector of unknown parameters. The vector of dependent variables (y), subject to general nonlinear, two-point boundary conditions (5) is approximated by a polynomial function (S(x)) on each subinterval [x_n, x_n+1] of a mesh a = x_o < x₁ < ··· < x_N=b.

This approximation should satisfy the boundary conditions (6) and satisfies the differential equations at both ends, as well at the midpoints of each subinterval (7) (8) (9)

In the context of this work, the integral term found on the right-hand side of Equation (1) was substituted by Gauss-Legendre quadratures 19. The collocation method is formally derived by evaluating the governing integro-differential equation at the collocation points, which results in a system of nonlinear ordinary differential–algebraic equations describing the evolution of the solution at the collocation points. This methodology is very attractive due to its easiness of implementation, even when the problem to be solved is highly nonlinear 19,20.

3. Mathematical formulation and solution of the inverse problem

In the inverse problem, we consider that the space-dependent albedo is unknown. The inverse problem would then consist on the estimation of the coefficients of this function. In the inverse formulation considered in this work the unknown space-dependent albedo is determined by using a parameter estimation approach. The albedo is thus estimated as a sampled function with a total of N_un discrete values and, therefore, we can write the vector of the unknowns as (10)

It is also considered that experimental data on the output radiation intensity measured at both boundaries, τ=0 and τ=τ_o, at different polar angles is available, i.e., Y_i, i=1, 2, … , N_d, where N_d is the total number of experimental data. From this experimental data, we then try to solve the inverse radiative transfer problem of estimating the space-dependent albedo.

As the number of experimental data, N_d, is assumed to be higher than the number of unknowns to be determined, N_un, the inverse problem is formulated as a finite dimensional optimization problem in which we seek to minimize the squared residues cost function (11) where I_i and Y_i are the calculated and measured values of the radiation intensity obtained at the same boundary, τ=0 and τ=τ_o, and at the same polar angle represented by μ_i, with i=1, 2, … , N_d. The elements of the vector of residues are given by (12)

As in any experimental procedure, there are uncertainties associated to the measured values. This aspect will be taken into account in this work by adding noise to the calculated solution. The radiative properties are known in this case. The standard deviation of experimental errors, σ, is included as shown in the following equation (13) where π_i is a pseudo-random number, generated with a zero mean Gaussian distribution and standard deviation equal to one, stands for the calculated values for the radiation intensity obtained by using the exact values for the unknowns.

4. DE algorithm

The DE algorithm is an optimization technique that belongs to the family of evolutionary computation and was proposed by Storn and Price 21,22. It is presented as a variant of genetic algorithms to achieve the goals of robustness in optimization and faster convergence to a given problem. DE differs from other evolutionary algorithms in the mutation and recombination phase and presents the following advantages: simple structure, easiness of use, speed, robustness and ability to escape from local optima 21,22.

DE creates new candidate solutions by combining the parent individual and several other individuals of the same population. A candidate replaces the parent only if it has a better fitness value. DE depends on three different parameters: the amplification factor of the difference vector (F), the crossover control parameter (CR) and the population size (NP).

The population of the canonical DE algorithm contains NP × D-dimensional vectors: (14) where G denotes the generation. During one generation for each vector, DE employs the mutation and crossover operations to produce a trial vector: (15)

Then, a selection operation is used to choose vectors for the next generation (G + 1). The initial population is selected uniformly randomly between the lower (x_L) and upper (x_U) bounds defined for each variable x_j. These bounds are specified by the user according to the nature of the problem. After the initialization, DE performs several vector transforms (operations), in a process called evolution.

4.1. Mutation operation

The mutation for each population vector creates a mutant vector: (16)

A mutant vector can be created by using one of the mutation strategies. The most useful strategies are the following:

1.	‘rand/1’: (17)
2.	‘best/1’: (18)
3.	‘current to best/1’: (19)
4.	‘best/2’: (20)
5.	‘rand/2’: (21)

where the indexes r₁, r₂, r₃, r₄ and r₅ represent the random and mutually different integers generated within the range [1, NP] and also different from index i. F is a mutation scale factor within the range [0, 2], usually less than 1. is the best vector in the generation G.

4.2. Crossover operation

After mutation, a ‘binary’ crossover operation forms the final trial vector, according to the i-th population vector and its corresponding mutant vector. (22) CR is a crossover parameter or factor within the range [0, 1] and presents the probability of creating parameters for trial vector from a mutant vector. Index j_rand is a randomly chosen integer within the range [1, NP]. It is responsible for the trial vector containing at least one parameter from the mutant vector. If the control parameters from the trial vector are out of bounds, the proposed solutions described in 21,22 are: reflected back within the bounds, set on the bounds or use as they are (out of the bounds).

4.3. Selection operation

The selection operation selects, according to the fitness value of the population vector and its corresponding trial vector, which vector will survive as a member of the next generation. For example, if we have a minimization problem, the following selection rule is used: (23) DE has been successfully applied to various fields such as digital filter design 23, estimation of heat transfer parameters in a bed reactor 24, synthesis and optimization of heat-integrated distillation systems 25, parameter estimation in fed-batch fermentation process 26, optimization of thermal cracker operation 27, engineering system design 15, apparent thermal diffusivity estimation during the drying of fruits 28, estimation of the parameters of Page's equation and heat loss coefficient by using experimental data from a realistic rotary dryer 29, treatment of multi-response surface applied to machinability of stainless steel using multi-objective optimization 30, determination of thermal properties in inverse heat conduction problems 31, turning process optimization 32, and other applications 22.

5. Results and discussion

In order to deal with more difficult test cases, we consider that only the left boundary of the medium at τ = 0 is subjected to the incidence of isotropic radiation originated at an external source while there is no radiation coming into the medium through the boundary at τ = τ_o. Therefore, the boundary conditions given by Equations (2) and (3) become: (24)

As real experimental data were not available, we generated synthetic experimental data by adding noise to the values calculated for the exit radiation intensities using the exact values of the radiative properties for each test cases (see Equation (13)). In all test cases, we have considered noiseless data as well as data with noise in the order of, or smaller than 10%.

In order to evaluate the performance of the DE algorithm we chose four representative test cases, with a space-dependent albedo given as piecewise constant or a polynomial function. The inverse solver estimates the space-dependent albedo as a set of 12 points, i.e. N_un equal to 12: (25)

For each set of experimental data (noiseless data and 5% and 10% error noisy data) the DE was executed using 20 different seeds for the random generation of the potential candidates. The results are presented always considering the average of the reconstructions for the 20 different seeds.

It should be emphasized that 20 points were used for the approximation of the variable μ (−0.9983, −0.9830, −0.9426, −0.8765, −0.7864, −0.6750, −0.5451, −0.4003, −0.2446, −0.0823, 0.0823, 0.2446, 0.4003, 0.5451, 0.6750, 0.7864, 0.8765, 0.9426, 0.9830, and 0.9983) and 10 collocation points (equally spaced) were taken into account to solve the direct problem. In all the test cases, internal detectors were considered for avoiding the non-uniqueness of the solution of the inverse problem.

5.1. Test case 1

In this case, the following function for the space-dependent scattering albedo is considered: (26)

Considering the function specified above, the values for the exit radiation intensities (I(τ, μ)) were calculated, totalizing 20 points for each approximation of the variable μ considered.

Table 1 presents the results obtained by DE algorithm considering various population sizes and 100 generations.

Table 1. Comparison of the DE performance considering various population sizes.

Noise level (%)	Population size (NP)	Objective function (Equation (11))	Neval
	10	1.999E-09	1010
0	25	1.222E-11	2525
	50	2.345E-12	5050
	10	4.099E-02	1010
5	25	2.223E-02	2525
	50	9.898E-03	5050
	10	0.394E + 01	1010
10	25	0.280E + 01	2525
	50	0.108E + 01	5050

In this table, as expected, the increase of population size (NP) led to the decrease of the objective function value, and consequently, the increase of the number of evaluations (N_eval) of the objective function.

Figure 2 presents the comparisons of the exact and estimated albedo using the DE algorithm with standard deviation of the experimental data equal to 0%, 5% and 10%.

Figure 2. Comparison of the exact and estimated space-dependent albedo using the DE algorithm considering 0%, 5% and 10% error in experimental data.

5.2. Test case 2

In this case, the following function for the space-dependent scattering albedo is considered: (27)

Considering the function specified above, the values for the exit radiation intensities were calculated (20 points for each approximation of the variable μ).

Table 2 presents the results obtained by DE algorithm considering various population sizes and 100 generations.

Table 2. Comparison of the DE performance considering various population sizes.

Noise level (%)	Population size (NP)	Objective function (Equation (11))	Neval
	10	1.089E-08	1010
0	25	3.011E-10	2525
	50	8.096E-11	5050
	10	1.039E-02	1010
5	25	1.039E-02	2525
	50	1.038E-02	5050
	10	0.422E + 01	1010
10	25	0.292E + 01	2525
	50	0.122E + 01	5050

In this table, also is observed that the increase of population size (NP) led to the decrease of the objective function value, and consequently, the increase of the number of evaluations (N_eval) of the objective function.

Figure 3 presents the comparisons of the exact and estimated albedo using the DE algorithm with standard deviation of the experimental data equal to 0%, 5% and 10%.

Figure 3. Comparison of the exact and estimated space dependent albedo using the DE algorithm considering 0%, 5% and 10%F error in experimental data.

5.3. Test case 3

In this case, the following function for the space-dependent scattering albedo is considered: (28)

As mentioned earlier, for the linear function defined by phases specified above, the values for the exit radiation intensities were calculated (20 points for each approximation of the variable μ).

Table 3 presents the results obtained by DE algorithm considering various population sizes and 100 generations.

Table 3. Comparison of the DE performance considering various population sizes.

Noise level (%)	Population size (NP)	Objective function (Equation (11))	Neval
	10	1.343E-10	1010
0	25	3.313E-11	2525
	50	2.222E-12	5050
	10	1.232E-02	1010
5	25	1.333E-02	2525
	50	9.878E-03	5050
	10	0.456E + 01	1010
10	25	0.343E + 01	2525
	50	0.129E + 01	5050

Figure 4 presents the comparisons of the exact and estimated albedo using the DE algorithm with standard deviation of the experimental data equal to 0%, 5% and 10%.

Figure 4. Comparison of the exact and estimated space dependent albedo using the DE algorithm considering 0%, 5% and 10% error in experimental data.

5.4. Test case 4

The next test case was proposed and solved by Stephany et al. 17 using the ant colony optimization (ACO) algorithm or its hybridization with the gradient-based Levenberg–Marquardt method. In this case, the following function for space dependent scattering albedo is considered: (29)

For the quadratic function specified, the values for the exit radiation intensities were calculated (20 points for each approximation of the variable μ).

Table 4 presents the value of the objective function obtained using three different strategies: pure ACO, with 120 ants and 500 generations; the hybridization between ACO and Levenberg–Marquardt, ACO–LM, with 10 ants and 40 generations and DE algorithm, with 10 individuals and 100 generations.

Table 4. Comparison of the performance for the three optimization schemes: pure ACO, hybridization ACO–LM and DE.

	Noise level (%)
Algorithm	0	5	10	Neval
ACO	7.22E − 06	–	–	60,000
ACO–LM	4.39E − 10	1.31E − 03	–	400 + 19
DE	1.638E − 09	1.272E − 05	2.343E − 01	1010

Comparing the evolutionary techniques, we can see that the DE algorithm achieves a better result with less computational effort when compared to the ACO algorithm and results of the same order of magnitude from those obtained by the ACO–LM method.

Figure 5 presents the comparison of the exact and estimated albedo using the DE algorithm with standard deviation of the experimental data equal to 0%, 5% and 10%, respectively.

Figure 5. Comparison of the exact and estimated albedo using the DE algorithm considering 0%, 5% and 10% error in experimental data.

In general, for the four test cases, we can say that the DE algorithm was able to obtain good results for the problem either with or without noise in experimental data.

6. Conclusions

In the present work, the DE algorithm was used for estimating space-dependent single scattering albedo in radiative transfer problems.

In this sense, four test cases were studied and it was possible to conclude that the methodology proposed in this work led to good results for an acceptable number of generations. It should be pointed out that the DE algorithm led to optimal values that are very similar to those obtained with the ACO algorithm or its hybridization with the gradient-based Levenberg–Marquardt method, requiring however a smaller number of objective function evaluations.

Another interesting aspect is that by adding noise to the synthetic experimental data results, as expected, an increase in the objective function values, as observed in all tables. Such, behaviour is expected because noise does not permit the convergence of the optimization evolution to the real experimental values. Consequently, the user should be aware of this behaviour when using real experimental data, which is always affected by noise.

As future research work, the authors intend to analyze the influence of changing the parameter values of the optimization algorithms and observe their influence on the performance of the estimation procedure.

Acknowledgements

Dr Lobato acknowledges the financial support provided by FAPEMIG, Fundação de Amparo à Pesquisa de Minas Gerais. Dr Steffen is thankful to the financial support provided by the agencies FAPEMIG and CNPq (INCT-EIE). Dr Silva Neto acknowledges the financial support provided by the research agencies FAPERJ, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro and CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico. This paper was presented at the Inverse Problems, Design and Optimization Symposium, João Pessoa, Brazil, 25–27 August, 2010.