Application of the generalized extremal optimization algorithm to an inverse radiative transfer problem

Abstract

The recently developed generalized extremal optimization (GEO) algorithm is applied for the solution of an inverse problem of radiative properties estimation. A comparison with two other stochastic methods, simulated annealing and genetic algorithms, is also performed, demonstrating that GEO is competitive. From the test case results we could also infer that a hybridization of GEO with gradient-based methods is very promising.

Keywords:

• Inverse problems
• Radiative transfer
• Optimization
• Generalized extremal optimization

1. Introduction

Several explicit and implicit formulations have been developed for the solution of inverse radiative transfer problems 1–3,, 8, 17, 19. In the implicit formulation, the inverse problem is usually replaced by an optimization problem in which we seek to minimize a cost function.

For the solution of the resulting optimization problem, deterministic, stochastic and hybrid approaches have been used 18, 20, 21. A growing interest has been observed towards heuristic methods because deterministic methods may not converge at all when a poor initial guess is used for the unknown parameters, and when they converge they may in fact lead to local minima 18. Recently, a new optimization algorithm inspired by a simplified evolutionary model was proposed, the so called generalized extremal optimization (GEO), 23. This algorithm was developed to be easily applicable to a broad class of nonlinear constrained optimization problems, with the presence of any combination of continuous, discrete and integer values, while having only one free parameter to be adjusted. Its efficacy to tackle complex design spaces has been demonstrated with test functions and real design problems 10, 24–27. Nonetheless, being a new algorithm, many of its features remain to be explored, such as parallelization, hybridization with other optimization algorithms, or different types of representation for the design variables. Moreover, it is also interesting to verify its efficiency in dealing with real problems coming from different areas of application, in comparison to other stochastic optimization methods.

In the present work, we apply the GEO for the solution of an inverse problem of radiative properties estimation. From the experimental data on the intensity of the exit radiation, we want to obtain estimates for the optical thickness, single scattering albedo, and the boundary diffuse reflectivities of one-dimensional homogeneous participating media. Test case results are presented, and comparisons with two other stochastic methods, simulated annealing (SA) and genetic algorithms (GA), are performed.

2. Mathematical formulation of the direct and inverse problems

2.1. Direct problem

Consider the problem of radiative transfer in an absorbing, isotropically scattering, plane-parallel grey medium with diffusely reflecting boundary surfaces. The mathematical formulation of the steady-state direct problem with azymuthal symmetry is given by 16: (1) (2) (3) where I(τ, μ) represents the dimensionless radiation intensity, τ the optical variable, μ the cosine of the polar angle, ω the single scattering albedo, and ρ₁ and ρ₂ the diffuse reflectivities at τ = 0 and τ = τ₀, respectively. The intensities of the external isotropic radiation sources are represented by A₁ and A₂.

When the geometry, the radiative properties and the boundary conditions are known, problem (1) may be solved yielding the values of the radiation intensity I(τ, μ), for 0 ≤ τ ≤ τ₀ and −1 ≤ μ ≤ 1. This is the direct problem. For the solution of the direct problem in the present work we use Chandrasekhar's discrete ordinates method 9.

2.2. Inverse problem

In the inverse radiative transfer problem we are interested in the following radiative properties which are considered unknown (4) but measured data on the intensity of the exit radiation at the boundaries τ = 0 and τ = τ₀, i.e., Y_i, i = 1, 2, … , N, are considered available, where N represents the total number of experimental data.

As the number of measured data, N, is usually much larger than the number of parameters to be estimated, the inverse problem is formulated as a finite dimensional optimization problem in which we seek to minimize the cost function (also referred to as the objective function) (5) where I_i represents the calculated value of the radiation intensity (using estimates for the unknown radiative properties ) at the same boundary, and at the same polar angle, for which the experimental value Y_i is obtained.

3. The generalized extremal optimization algorithm (GEO)

The GEO algorithm 24–26, is a newly proposed meta-heuristic approach suitable to tackle complex optimization problems. It was developed as a generalization of the extremal optimization (EO) method 6, in such a way that it can be readily applied to a broad class of problems. Either EO or GEO are based on the simplified evolutionary model of Bak–Sneppen 5, which was devised to show the emergence of self-organized criticality (SOC) in ecosystems. The theory of SOC has been used to explain the power law signatures that emerge from many complex systems in such different areas as geology, economics, and biology 4. It states that large interactive systems evolve naturally to a critical state where a single change in one of its elements generates “avalanches” that can reach any number of elements in the system. The probability distribution of the sizes “s” of these avalanches is described by a power law in the form (6) where γ is a positive parameter. That is, smaller avalanches are more likely to occur than big ones, but even avalanches as big as the whole system may occur with a nonnegligible probability. This kind of dynamic behaviour, according to Bak and Sneppen 5, could explain the bursts of evolutionary activity observed in the fossil record and that has been given the name of punctuated equilibrium 12.

In the Bak–Sneppen model, species are represented in a lattice and, for each of them, there is associated a fitness number in the range [0, 1]. The evolution is simulated forcing the least adapted species, the one with the least fitness, and their neighbours, to change (it can “evolve” or be “extinct” and replaced by a new one, that not necessarily has a better fitness). This is done by assigning new fitness numbers, randomly, to these species. This simulated ecosystem is started with the fitness of the species distributed uniformly in the range [0, 1]. Since the less adapted species are constantly forced to change, the average fitness value of the ecosystem increases and, eventually, some time after initialization all species have a fitness number above a “critical level”. However, as even good species may be forced to change (if they are neighbours of the least adapted one), it happens that a number of species may fall below the critical level from time to time. That is, the equilibrium (being above the critical level in “stasis”) of one or more species is punctuated by avalanches, whose occurrence is described by a power law 5. Although the claim that the evolution of species may happen in a system that is self organized critical has been controversial 13, 22, an optimization heuristic based on the Bak–Sneppen model may evolve solutions quickly, systematically mutating the worst individuals, preserving at the same time throughout the search process, the possibility of probing different regions of the design space (via avalanches).

Based on the Bak–Sneppen model, Boettcher and Percus developed the EO method 6, which has been applied successfully to combinatorial optimization problems 6, 7, in which a fitness number is associated with the design variables. However, as they pointed out, in some cases this may become an ambiguous or even impossible task 6. The GEO algorithm was devised to overcome this problem, so that it could be easily applicable to a broad class of problems, with any kind of design variable, either continuous, discrete or a combination of them, in a design space that may be multimodal or even discontinuous and subject to any kind of constraint.

In GEO a string of L bits is considered a population of species. That is, each bit is a species. The string encodes the M design variables. For each of them is associated a fitness number that is proportional to the gain (or loss) the objective function value has in mutating (flipping) the bit. All bits are then ranked from 1, for the least adapted bit, to L for the best adapted. A bit is then mutated according to the probability distribution P ∝ k^−γ, where k is the rank of a selected bit candidate to mutate, and γ is a free control parameter. Making γ → 0, the algorithm becomes a random walk, while for γ → ∞, we have a deterministic search. It has been observed that the best value of γ, i.e., the one that yields the best performance of the algorithm, for a given application generally lies within the range [0.75, 3.0]. This, added to GEO having only one free parameter, makes the algorithm easily settable to give its best performance on the problem that is being tackled. After the bit is mutated, the procedure is repeated until a given stopping criterion is reached, and the best configuration of bits (the one that gives the best value for the objective function) found is returned.

In a variation of the canonical GEO described earlier, the bits are ranked separately for each sub-string that encodes each design variable, and M bits (one for each variable) are flipped at each iteration of the algorithm. In previous works, 24, 25, it was observed that this implementation seems to be more efficient than the canonical one for the cases in which the problem being tackled has only bound constraints (constraints that represent the limits of feasible values for the design variables). Since this is the case here, this approach was implemented to tackle the inverse radiative transfer problem described in section 2. The main steps of the implementation are shown in figure 1.

Figure 1. Outline of the GEO algorithm as implemented for the solution of the inverse radiative transfer problem.

The performance of GEO in each particular application of interest is affected by the value of its free parameter γ. Therefore, a preliminary study was done so that a proper value of γ could be used in the inverse radiative transfer problem. For each value of γ taken from a discrete set [0.75, 2.25] GEO was run 20 times, for a test case considering noiseless data, with the following radiative properties: τ₀ = 1.0, ω = 0.7, ρ₁ = 0.1 and ρ₂ = 0.1. For computationally expensive objective functions, such as the one investigated here, the search for the proper value of γ is usually done with fewer function evaluations per run than the value intended to be used later on the solution of the problem at hand. For the preliminary study each run was stopped after 10,000 function evaluations. The average of the best values of the objective function found in the runs for each γ are shown in figure 2.

Figure 2. Average of the best objective function values found in 20 runs of GEO for each value of γ.

From figure 2 it can be seen that the best performance of GEO for this problem occurs for γ close to 1.25. This value was used for all other runs performed in the present study, since it is not expected that the best value of parameter γ, i.e., the one that yields the best performance for a given application, changes significantly for different instances of the same problem.

A detailed explanation of GEO can be found in 25, 26, while examples of applications to real complex design problems can be found in 10, 25–27.

4. The simulated annealing and genetic algorithms

In order to evaluate the performance of GEO, we present in this work a comparison with the results obtained with two other stochastic methods, SA and GAs. We present next just a brief description of such methods. Comprehensive descriptions or further details of the computational implementation may be found in 11, 14, 1518, 20.

4.1. Simulated annealing

Based on statistical mechanics reasoning applied to a solidification problem Metropolis et al. 15 introduced a simple algorithm that can be used to accomplish an efficient simulation of a system of atoms in equilibrium at a given temperature (T). In each step of the algorithm a small random displacement of an atom is performed and the variation of the energy ΔE is calculated. If ΔE < 0 the displacement is accepted, and the configuration with the displaced atom is used as the starting point for the next step. In the case of ΔE > 0, the new configuration can be accepted according to Boltzmann probability (7) where k_B is the Boltzmann's constant.

A uniformly distributed random number r in the interval [0, 1] is then calculated and compared with P(ΔE ). The Metropolis criterion establishes that the new configuration is accepted if r < P(ΔE), otherwise it is rejected and the previous configuration is used again as a starting point. Using the objective function , defined in equation (3) in place of energy, and defining configurations by a set of variables {Z_i}, i = 1, 2, … , M, see equation (2) where M = 4 is the total number of unknowns, the Metropolis procedure generates a collection of configurations of a given optimization problem at some temperature T 14. This temperature is simply a control parameter. The simulated annealing process consists of first “melting” the system being optimized at a high effective “temperature”, then lowering the “temperature” until the system “freezes” and no further change occurs.

The main control parameters of the algorithm implemented (“cooling procedure”) are the initial “temperature”, T₀, the cooling rate, r_t, the number of steps performed through all elements of vector , N_S, the number of times the procedure is repeated before the “temperature” is reduced, N_t, and the number of points of minimum (one for each temperature) that are compared and used as the stopping criterion if they all agree within a tolerance ϵ, N_ϵ.

4.2. Genetic algorithms

Darwin's theory of survival of the fittest gives the main idea of the method 11. A set of feasible designs constitutes a generation, which has a fixed number of individuals. A set of better designs is derived from the previous generation where the individuals are allowed to reproduce and cross among themselves with bias allocated to the fittest members. Combinations of the most favourable characteristics of the mating members of the population result in a new generation that is more fit than the previous one.

In the traditional GA each set of variables is represented by a binary string, and the cost function , see equation (3), is used to construct the fitness function that indicates how good a member is in his generation.

The GA is implemented with three basic operations: reproduction, crossover and mutation, and the main control parameters of the algorithm implemented are the number of individuals in the population, n, probability of crossover, p_c, probability of mutation, p_m, and number of generations, n_g. In this study we used a micro-genetic algorithm (μGA), which is a variation of the simple GA. The major difference between the traditional GA and the μGA comes in the population choice and the way new information is brought into the evolution process. Usually the population size is fixed at five. It is known that GA generally do poor with very small populations due to insufficient information processing and early convergence to nonoptimal solutions. In the case of the μGA new strings are brought into the population at regular intervals. The restart procedure helps in avoiding premature convergence.

5. Results and discussion

In order to evaluate the performance of the algorithms GEO, SA, and μGA in the solution of the inverse radiative transfer problem described in section 2, the three sets of radiative properties shown in table 1 were considered. These sets were chosen for yielding relatively difficult test cases for the evaluation of the algorithms. The following intensities for the external radiation sources were considered: A₁ = 1.0 and A₂ = 0.0.

Table 1. Exact values of the radiative properties

Radiative property	Case 1	Case 2	Case 3
Optical thickness, τ0	1.0	2.0	0.5
Single scattering albedo, ω	0.5	0.8	0.3
Diffuse reflectivity, ρ1	0.2	0.1	0.1
Diffuse reflectivity, ρ2	0.2	0.8	0.8

As real experimental data on the intensity of the exit radiation was not available, we have generated sets of synthetic experimental data with (8) where represents the calculated values of the radiation intensity using the exact values of the radiative properties, , as given in table 1, r_i is a pseudo-random number in the range [−1, 1], σ simulates the SD of the measurement errors, and as described in section 2.2, N represents the total number of experimental data. The values of σ = 0.005, 0.002 and 0.0025 lead to errors in the order of, or smaller than, 5% in the exit radiation intensives for Cases 1, 2 and 3, respectively.

As mentioned before, in section 2.1, Chandrasekhar's discrete ordinates method is used for the solution of the direct radiative transfer problem. A Gauss–Legendre quadrature is then used to replace the integral term in equation (1a). Therefore, the collocation points correspond to the roots of the Legendre polynomials. These values are then related to the cosine of the polar angles, μ_n, with n = 1, … , N, being half of them, n = 1, … , N/2 in the positive range, μ > 0, and half of them, n = N/2 + 1, … , N, in the negative range, μ < 0. Half of the experimental data is acquired at τ = τ₀, i.e., Y_i, with i = 1, … , N/2, at the same polar angles for which μ_n corresponds to μ > 0, i.e n = 1, … , N/2, and half of the experimental data is acquired at τ = 0, i.e., Y_i, with i = N/2 + 1, … , N, at the same polar angles for which μ_n corresponds to μ < 0, i.e., n = N/2 + 1, … , N.

In figure 3 is presented the evolution of the average of the best values of the objective function, in 10 runs, and for each method (GEO, SA and μGA), for Case 1 listed in table 1, using experimental data without noise, i.e., σ = 0 in equation (6). In figure 4 the same test case is considered, but now with noisy data, σ = 0.005.

Figure 3. Average of the best values of the objective function, as a function of the number of function evaluations for Case 1, without noise.

Figure 4. Average of the best values of the objective function, as a function of the number of function evaluations for Case 1, with noise.

In table 2 are presented for Case 1, with and without noise in the synthetic experimental data, the worst, average and best estimates obtained for the unknown radiative properties (design variables). Here the worst estimates obtained for each method correspond to the run, among the 10 runs performed for each method, for which the objective function is the highest at the end of the run, and the best estimates correspond to the run for which the value of the function is the lowest.

Table 2. Worst, average and best estimates for Case 1

			τ0	ω	ρ1	ρ2	fx^a
Exact value			1.000	0.500	0.200	0.200	–
SA	Without noise	Worst	1.000	0.499	0.198	0.200	6.07 × 10^−8
		Average^b	1.000	0.500	0.199	0.200	1.80 × 10^−8
		Best	1.001	0.500	0.201	0.200	1.84 × 10^−9
	With noise	Worst	1.032	0.525	0.253	0.204	4.03 × 10^−4
		Average^2	1.034	0.526	0.254	0.203	4.03 × 10^−4
		Best	1.034	0.526	0.254	0.203	4.03 × 10^−4
μGA	Without noise	Worst	0.963	0.509	0.251	0.250	2.50 × 10^−4
		Average	1.010	0.508	0.218	0.201	4.93 × 10^−5
		Best	0.987	0.495	0.190	0.205	2.79 × 10^−6
	With noise	Worst	1.540	0.817	0.000	0.000	2.2755
		Average	1.505	0.780	0.092	0.019	1.8205
		Best	1.173	0.589	0.373	0.163	5.58 × 10^−4
GEO	Without noise	Worst	0.943	0.497	0.231	0.253	2.63 × 10^−4
		Average	0.965	0.481	0.157	0.211	1.31 × 10^−4
		Best	1.012	0.502	0.200	0.194	9.00 × 10^−6
	With noise	Worst	0.924	0.503	0.234	0.266	9.46 × 10^−4
		Average	1.030	0.527	0.256	0.212	6.16 × 10^−4
		Best	1.046	0.541	0.295	0.212	4.37 × 10^−4
^afx is the value of the objective function at the end of the run.
^bThe SA stopping criterion was reached for all runs in <25,000 evaluations.

In figures 5 and 6, and in table 3, the results for the Case 2 listed in table 1 are presented, and in figures 7 and 8, and table 4, the results for Case 3 are presented.

Figure 5. Average of the best values of the objective function, as a function of the number of function evaluations for Case 2, without noise.

Figure 6. Average of the best values of the objective function, as a function of the number of function evaluations for Case 2, with noise.

Figure 7. Average of the best values of the objective function, as a function of the number of function evaluations for Case 3, without noise.

Figure 8. Average of the best values of the objective function, as a function of the number of function evaluations for Case 3, with noise.

Table 3. Worst, average and best estimates for Case 2

			τ0	ω	ρ1	ρ2	fx^a
Exact value			2.000	0.800	0.100	0.800	–
SA	Without noise	Worst	2.003	0.800	0.101	0.800	1.72 × 10^−8
		Average^b	2.000	0.800	0.100	0.800	8.64 × 10^−9
		Best	2.000	0.801	0.102	0.800	0.92 × 10^−9
	With noise	Worst	2.027	0.803	0.109	0.799	7.76 × 10^−5
		Average^b	2.017	0.801	0.104	0.799	7.75 × 10^−5
		Best	2.016	0.801	0.104	0.799	7.75 × 10^−5
μGA	Without noise	Worst	2.845	0.851	0.243	0.732	1.66 × 10^−4
		Average	2.240	0.821	0.169	0.789	7.11 × 10^−5
		Best	1.955	0.790	0.062	0.797	7.63 × 10^−6
	With noise	Worst	3.847	0.894	0.402	0.669	4.01 × 10^−4
		Average	3.751	0.878	0.336	0.602	3.10 × 10^−4
		Best	2.933	0.853	0.249	0.720	2.50 × 10^−4
GEO	Without noise	Worst	4.844	0.891	0.373	0.249	4.28 × 10^−4
		Average	2.507	0.817	0.148	0.718	1.72 × 10^−4
		Best	2.033	0.809	0.138	0.803	1.80 × 10^−5
	With noise	Worst	2.077	0.816	0.164	0.831	5.21 × 10^−4
		Average	2.385	0.824	0.167	0.774	2.62 × 10^−4
		Best	2.160	0.811	0.129	0.787	9.90 × 10^−5
^afx is the value of the objective function at the end of the run.
^bThe SA stopping criterion was reached for all runs in <25,000 evaluations.

Table 4. Worst, average and best estimates for Case 3

			τ0	ω	ρ1	ρ2	fx^a
Exact value			0.500	0.300	0.100	0.800	–
SA	Without noise	Worst	0.500	0.300	0.100	0.800	3.08 × 10^−8
		Average^b	0.500	0.300	0.100	0.800	1.57 × 10^−8
		Best	0.500	0.300	0.100	0.800	8.49 × 10^−10
	With noise	Worst	0.487	0.305	0.147	0.806	1.08 × 10^−4
		Average^b	0.486	0.305	0.149	0.806	1.08 × 10^−4
		Best	0.487	0.305	0.148	0.806	1.08 × 10^−4
μGA	Without noise	Worst	0.430	0.362	0.375	0.841	5.82 × 10^−4
		Average	0.463	0.319	0.226	0.818	2.08 × 10^−4
		Best	0.518	0.299	0.049	0.793	2.52 × 10^−5
	With noise	Worst	0.254	0.343	0.680	0.902	3.22 × 10^−3
		Average	0.431	0.327	0.298	0.833	8.01 × 10^−4
		Best	0.469	0.325	0.233	0.818	1.6 × 10^−4
GEO	Without noise	Worst	0.396	0.210	0.152	0.750	1.91 × 10^−2
		Average	0.501	0.290	0.107	0.786	4.32 × 10^−3
		Best	0.523	0.277	0.003	0.789	2.08 × 10^−4
	With noise	Worst	0.831	0.543	0.094	0.728	3.31 × 10^−2
		Average	0.525	0.378	0.237	0.812	6.82 × 10^−3
		Best	0.386	0.302	0.382	0.844	8.17 × 10^−4
^afx is the value of the objective function at the end of the run.
^bThe SA stopping criterion was reached for all runs in <25,000 evaluations.

From the results presented in figures 3–8, it can be seen that the SA performed on average better than the μGA and GEO for all cases, either with or without noise. Nonetheless, it is interesting to note that the SA usually did not have the best performance early in the search. In fact, its results became the best only after around 7500 function evaluations. This is probably the consequence of the annealing, and it is an indication that, for this problem, making the algorithm more deterministic along the search is beneficial. It also forecasts that the performance of GEO could be significantly improved if the parameter γ is varied through the search, which would have an effect similar to the decrease of temperature in the SA. This will be a subject for a future study.

Another interesting result from figures 3–8 is that while GEO performed worse than μGA for the cases without noise, when noise was added it had a similar performance for Case 2 and a much better one for Case 1. In fact, for Case 1 with noise GEO performed close to the SA (even better when the number of function evaluations was <7000), while the performance of μGA was severely degraded. The μGA had a better performance than GEO for the instances with noise only in Case 3, the one with the lowest level of noise. That is, GEO seems to be more robust to noise than μGA.

As described before, 10 runs were performed for each case (1–3 in table 1), for each method (GEO, SA and μGA), considering two sets of experimental data for each case, i.e., with and without noise. Each run for GEO and μGA was initiated from a randomly chosen initial guess. Both GEO and μGA used 10 bits to encode each of the unknowns, which implies a precision of 0.005 for τ, and 0.001 for ω, ρ₁ and ρ₂. For both methods the computations were stopped at 25,000 function evaluations. The choice of the control parameter for the GEO, γ, has already been described.

The control parameters used for the SA were: T₀ = 5.0, r_t = 0.75, N_s = 20, N_t = 5, N_ϵ = 4, and the convergence tolerance was ϵ = 10⁻⁶.

The runs of GEO were performed on one PC with the processor Intel Xeon (2.2 GHz with 1 GB of RAM) and one PC with the processor AMD Athlon (1.8 GHz with 256 MB of RAM). Each run of GEO took an average of 2 h 48 min on the former processor and 2 h 38 min on the latter. The runs of SA and μGA were performed on one PC with the processor Pentium III (750 MHz with 192 MB of RAM). On average each SA run took 3 h 25 min and each μGA run took 3 h 10 min.

The better values obtained for the objective function with SA are due to the estimated values for the radiative properties, which are closer to the exact values than the ones obtained with the μGA and GEO, as shown in tables 2–4. However, even for the SA, the presence of noise can degrade the performance of the algorithms significantly, mainly on the recovery of the value of the parameter ρ₁. This is in fact the main difficulty associated with the test cases presented, i.e., the low values of ρ₁ combined with the external illumination given by A₁ = 1.0 and A₂ = 0.0, in equation (1b) and (c), respectively. The external radiation detectors receive the information regarding the diffuse reflectivity at the inner part of the boundary at τ = 0 only after, the radiation goes into the medium through this boundary, is reflected at the boundary τ = τ₀. Then it is partially transmitted through the boundary τ = 0, being then measured by external detectors located on this boundary, and then partially reflected, being then partially transmitted through the boundary τ = τ₀, when it is measured by the external detectors located on that side of the medium. The low sensitivity of the exit radiation intensities to the diffuse reflectivity ρ₁ is confirmed by a sensitivity analysis.

As shown in 18, 20, a hybrid method which combines a global search metaheuristic and a gradient-based method (GBM), would most likely yield better results for the inverse radiative transfer we are looking at. It is interesting to note that in a hybrid strategy where the metaheuristic is used for a few number of function evaluations, just to find regions of local optima to be exploited by the GBM, the use of GEO or GA could yield a better combination than with the SA since, in most cases, they are more efficient than the SA in reducing the value of the objective function at the beginning of the search.

6. Conclusions

In the present work, it is demonstrated that the recently developed GEO algorithm performed well, in a competitive basis with other stochastic methods, namely SA and GA, for the solution of the inverse problem of radiative properties estimation.

It was devised also that a very promising approach would be the hybridization of GEO with a GBM. Another path to be followed in the near future is the use of a varying control parameter for GEO, γ, in a similar way to the cooling rate of SA.

The test cases considered in the present work were intentionally difficult, and the three methods used, GEO, SA and μGA performed competitively, with a slight advantage for GEO and μGA at the beginning of the search (what is an interesting feature for the hybridization with a GBM), and for SA in the long run (in case the stochastic algorithms are run fully).

Acknowledgments

The authors acknowledge the financial support provided by CNPq – Conselho Nacional de Desenvolvimento Científico e Tecnológico, CAPES – Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, and FAPERJ - Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro.