A specialized genetic algorithm for the electrical impedance tomography of two-phase flows

TECHNICAL PAPERS

A specialized genetic algorithm for the electrical impedance tomography of two-phase flows

Vanessa P. Rolnik^I; Paulo Seleghim Jr.^II

^[e-mail] Faculty of Philosophy, Sciences and Letters of Ribeirao Preto University of Sao Paulo - USP 14040-901 Ribeirão Preto, SP. Brazil

^IISenior Member, ABCM [e-mail] School of Engineering of Sao Carlos University of Sao Paulo - USP 400 13560-970 São Carlos, SP. Brazil

ABSTRACT

The main objective of this work is to contribute to the development of a new tomographic reconstruction method well suited for processing signals obtained from electrical or other soft sensing field probes. The adopted approach consists in formulating the reconstruction problem in terms of an error function, which assesses the difference between a prospective and the actual internal contrast distribution (3D image), and searching for its minimum with the help of a specialized genetic algorithm (GA). Numerical simulations have been performed to demonstrate the feasibility of the proposed reconstruction method, as well as to emphasize the relation between the ill-posed nature of the problem and the topology of the minimization hyper-surface, and the importance of considering this relation when designing the numerical solution procedure. Results show that convergence is greatly enhanced when a priori information is introduced in the error function.

Keywords: Electrical impedance tomography, genetic algorithm, inverse problem, optimization, multiphase-flow

Introduction

Industrial process tomography, and particularly applications involving multiphase flows, should be based on inexpensive and sufficiently robust sensing techniques to be applicable in uncontrolled environments. For instance, despite the quality of the images obtained with the use of X rays, positron emission and nuclear magnetic resonance, the use of such equipment in the harsh conditions of a fluidized coal combustor or an underwater petroleum pipeline is highly problematic. The natural choices usually fall within techniques based on electrical impedance measurements (conductive, capacitive or inductive) or on acoustical measurements (transmission, scattering, time-of-flight), mainly due to their low cost and robustness. The intrinsic problem related to these sensing techniques is that the acquired data are related to the sensed medium through strongly ill-posed differential/integral operators. In addition to this, it is known that reconstruction techniques based on linearization or some other artificial regularization procedure achieve stability at the cost of resolution and distinguishability, with possible spurious artifacts which need to be corrected a posteriori (Lányi, 1998; Knowles, 1998 and Borcea, 2003).

One of the main purposes of this work is to contribute to the development of new reconstruction algorithms well suited for soft field sensing techniques, and particularly for electrical impedance tomography (EIT). Our basic idea is to start from a prospective contrast distribution of the sensing volume, based for instance on the qualitative images delivered by a direct imaging probe (Seleghim and Hervieu, 2003), and to refine it iteratively until the predicted variables match the corresponding experimental measurements. In an EIT problem this can be done by devising an error function assessing the difference between the actual and the prospective electrical contrast distribution (conductance or permittivity), through the comparison of boundary measurements (electric current or charge distribution respectively) with the corresponding predictions from the model. The reconstruction procedure becomes thus a nonlinear minimization problem on the coefficients of a parameterized or a discretized version of the prospective internal contrast. Parameterization can be achieved according to many possible approaches such as, for instance, Fourier (Calderon, 1980; Allers and Santosa, 1991) or wavelet basis (Dobson, 1992). Another possible alternative is to restrict the contrast distribution on a discretization grid of the sensing volume. Although an adaptive grid would certainly be an optimal choice, as discussed by Cherkaeva and Tripp (1996), an equidistant grid has the characteristic of distributing errors more evenly.

Another problem that we intend to approach in this work is the design of an efficient and robust numerical minimization procedure for the nonlinear EIT problem. Even though existence and uniqueness of solutions hold for a large class of practical situations (Kohn and Vogelius, 1984, Kohn and Vogelius, 1985, Nachman, 1988 and Nachman, 1996), if no a priori restriction is imposed the inverse relation between excitation-response and the internal contrast distribution is discontinuous (Hadamard, 1902). Consequently, the problem becomes extremely sensitive and probably impossible to be solved within realistic experimental conditions, where boundary values are not fully known, measurements are noisy, calculations are subjected to round-off errors, etc. Numerical regularization techniques are thus mandatory in order to palliate the intrinsic ill-posed mathematical characteristic of an EIT problem. There is, in fact, extensive literature on such numerical procedures which assure convergence of reconstruction algorithms by enforcing some a priori information about the internal contrast, at the eventual cost of introducing spurious artifacts in the images (Engl et al., 1996). However, the question of which regularization methods are more effective regarding a given practical situation in multiphase flow EIT problems constitutes an open and interesting research area.

In addition to the high sensitivity to noise and round-off errors, the intrinsic ill-posed nature of the problem is also associated with irregularities on the optimization hyper-surfaces, the so-called pathology of the optimization problem, which can mislead or restraint the convergence of the numerical procedure. Among such pathological irregularities, which generally surround the global minimum, are multiple local minima (which is characteristic of void fraction probes inverse calibration problems, helical-shaped hyper-channels, hyper-saddles or hyper-plateaus (Rolnik and Seleghim, 2002). Convergence to the correct solution in such conditions is extremely difficult. For instance, the minimization procedure will most likely converge to the wrong minimum, in the case of multiple local minima, if initial guesses are not close enough to the global solution. Or else, if the correct solution is located inside a helical channel a marching type optimization procedure will converge only if extremely small refinement steps are executed, in which situation the overall computational time may become prohibitive. It is important to stress that both excitation/measurement strategies and the form of the error function, with or without regularization terms, have a major influence in these pathological features (Cherkaeva and Tripp, 1996; Figueroa and Seleghim, 2001).

Under these circumstances a reliable and efficient optimization procedure should be able to handle such problematic topological features, which are a consequence of the ill-posed nature of the inverse EIT problem, and, also, of converging rapidly to the correct solution if on-line application of the method is envisaged. A possible strategy to obtain such performance is to combine the qualities of a directional type method to the virtues of a random search method. The first is necessary to accelerate convergence in regular regions of the optimization hyper-surface, while the last is extremely important, for instance, to escape from local minima, to correct a wrong convergence trajectory or to progress on a flat region of the optimization hyper-surface. Important classes of optimization methods that comply with these performance requirements are the so-called evolutionary methods which, in a certain way, mimic the biological evolution process and Charles Darwin's natural selection theory (survival of the fittest) for solving complex and/or high dimensional optimization problems. In agreement with these ideas a specialized genetic algorithm (GA) is proposed in this work to solve the inverse nonlinear EIT problem. Among many interesting characteristics, which contrast with most classical optimization methods, a GA requires no local gradient information and, since the search for the minimum is performed within populations of possible solutions, it can be easily parallelized to handle multiple regularized variants of the error function. Furthermore, as it will be shown in next sections, a GA can incorporate in very a straightforward manner a priori information about the problem, such as bounds on the optimization parameters (local contrasts), phase fraction (global contrast), symmetry, etc. (Hsiao et al., 2001).

Statement of the Problem

The governing non-dimensional equations of an electrical impedance tomography problem can be derived from Maxwell equations according to the following:

where f represents the electric field, s the medium's contrast (permittivity, conductivity or permeability), f_exc (x, h) the excitation voltage profile and Q(x, h) the measured profile (electric charges or current). The variables f_exc and Q represent the stimulus and response respectively, and shall be used to determine f and to reconstruct the internal contrast s within the sensing volume W. It is also possible to adopt a Neumann sensing strategy, based on injecting currents and measuring voltage profiles along the boundary, but this usually produces data poorly correlated with the internal contrast distribution.

In mathematical terms, stimulus and response are canonically conjugated to the internal contrast through a formal operator (direct problem)

The reconstruction problem can thus be interpreted as inverting with respect to s, that is (inverse problem)

The dimensionality involved in the last expressions was deliberately made explicit to stress the following point: while in the direct problem (4) a 2D profile and a 3D distribution generate another 2D profile, in the inverse problem (5) two 2D profiles generate a 3D distribution. In other words, in the direct problem the information originally distributed over 5 dimensions is compressed into 2 dimensions with expected degradation or loss. Thus, the response Q(x, h) can be seen as an imperfect projection of the internal contrast s(x,y,z). In the inverse problem this imperfect Q(x,h) is used to reconstruct the s(x, y,z), and the problem is clearly ill-posed.

This being, let s_actual be the actual electric contrast within W to be reconstructed by refining a prospective (or an approximated) contrast distribution denoted s_approx. The corresponding responses to a common excitation profile f_exc (x, h) are denoted respectively Q_approx (x, h) and Q_actual (x, h), the first being determined through a numerical model of the problem and the later by direct measurement. According to (4) this can be expressed mathematically by the following

The difference between s_actual and s_approx (as seen through a particular perspective defined by f_exc) can be quantified through the following error function

In the above expressions s_actual is not known a priori and Q_actual must be determined by experimental measurements, while Q_approx can be formally calculated through (Eq. 6) since s_approx and f_exc are both known. More precisely, Q_approx is determined in two steps:

1 – calculate the electrical field f in W by solving a numerical version of Eq. 1 (finite difference, finite elements, etc.) with s = sapprox and f_exc (x, h) imposed on the boundary according to Eq. 2 and

2 – use the solution f and s_approx in expression (3) to calculate Q_approx (x, h). Once Q_actual and Q_approx have been determined, the error between them can be assessed through expression (8).

Within this framework, the reconstruction procedure consists simply in searching for the approximated contrast distribution s_approx that minimizes the error (8), so that a good match between Q_actual and Q_approx is achieved.

From a practical point of view, although theoretically a single pair of conjugated distributions f_exc and Q_actual is sufficient to determine the corresponding contrast distribution s_actual – theorem 5.1 in Nachman (1988) – the presence of noise and unavoidable experimental errors can produce inconsistent results depending on the sense that the problem is being solved. In the direct problem (4) the errors are strongly damped and the result is quite close to the correct one. On the other hand, in the inverse problem (5), the errors are greatly amplified and the reconstructed contrast distribution may be completely corrupted.

This issue has been studied in previous works. For instance, in the problem of determining bubble sizes from pierced length histograms, Seleghim Jr. and Milioli (2001) adopted a similar approach as the one adopted here and observed that relative errors of less than 0.001% in the input data could distort the calculations, producing astonishing deviations of more 1000% in the reconstructed distributions. A plausible strategy to handle this problem, already tested with success in ill-posed integral equations associated to inverse calibration problems (Seleghim Jr. and Schiavon, 2001), is to execute redundant measurements and to consider them altogether in a same error function. For errors and noise having a random ergodic behavior it is very likely that the information degraded in one particular test may be recovered through other experiments. Let thus {fⁱexc(x, h)}, for i=1,2...N, be a sequence of different excitation profiles sequentially applied on the boundary surface in a time interval that may be considered negligible compared with the evolution of the contrast distribution s due to the flow. To each fⁱexc(x, h) it is possible to associate two response profiles by direct measurements and by assuming a common s_approx. Synthetically this can be written as:

Hence, a global error function can be constructed from the individual errors

and be expressed by

if the amplitudes are to be emphasized or, alternatively,

if both amplitudes and inclinations should be simultaneously considered. These types of global error functions allow not only the solution of the inverse problem within realistic error conditions, but also contributes in regularizing some of the most negative features of the optimization surface, such as saddle points, sharp ridges, boundary minima, etc. These and other points will be illustrated by the numerical experiments presented in the sequel.

Numerical Simulations and Error Surface Analysis

To demonstrate the feasibility of the proposed reconstruction method let us consider the problem depicted in the Fig.1, in which the axes x y z describe the sensing volume W and (x,h) describes its boundary ¶W. A discrete version of the governing Eq. 1 can be derived in Cartesian coordinates, according to a central finite difference scheme:

where the coefficients are given by

The implicit notation is the following

The limits in (15) where fixed at N_x = 20, N_y= 20 and N_z= 60 after a convergence analysis.

Boundary conditions were adopted to simulate a common practical situation in which the longitudinal faces are grounded except in a single point where an excitation voltage of 10 V is applied (Dirac excitation). Null flux conditions are applied on the transversal sections to limit the extension of the sensing volume. The actual permittivity distribution was defined to represent the presence of a small 3D inclusion occupying a unitary mesh volume at the position indicated in Fig.1. Contrast values were chosen to model an air-water mixture, that is s_actual = 1 within the inclusion (air) and s_actual = 80 outside of it (water). To be able to exhibit the irregularities of the optimization hyper-surface a 2-dimensional problem was first studied. Different contrast approximations s_approx were generated by translating the prospective inclusion along the horizontal and vertical planes, indicated in gray in Fig.1. These contrast distributions were used in (4) (with the same excitation profile) to estimate Q_approx (x, h)½_{m, n} and n defining the position of the inclusion as indicated in Fig.1. Subsequently, error values e(m,n) were calculated with (8), in which Q_actual (x, h) was obtained with the same excitation profile and the actual contrast distribution. The results are shown in Fig.2.

We can clearly observe that these error surfaces exhibit an extremely pronounced minimum where the inclusions match each other, i.e. when s_actualºs_approx, indicating that the reconstruction can be accomplished by finding the global minimum of (8). This approach has several advantages. First, there are no restrictive hypotheses, such as the sensing field being 2-dimensional and parallel, what is not likely to occur in high contrast flows even with the use of confinement techniques (guard electrodes). In addition, this approach is virtually free of physical averaging effects and intrinsically 3D, provided the modeled regions upstream and downstream from the excitation and measurement electrodes are sufficiently long. Furthermore, the formulation of the problem is so that the introduction of additional information that contributes to the reconstruction process can be accomplished in a very straightforward manner. This last point, probably the most promising implication, concerns the development of the so called multimodal tomography based, for instance, on simultaneous electrical and acoustical measurement, in which, hypothetically, an object not seen by one would be captured be the other sensing field.

An analysis of Fig.2 also reveals the potential problems related to the topology of the error surface. In particular, the existence of saddle points, boundary minima and, above all, the nearly flat region surrounding the solution minimum (plateau), can seriously compromise the effectiveness of a numerical minimization methods. For example, a marching method based on the local inclination would inevitably stop if the optimization trajectory falls into this plateau. Otherwise, the same minimization procedure would converge to some point on the contour if the starting point had been set between the ridges and the external boundaries. This problem can be interpreted as the manifestation of the ill-posed nature of the inverse problem as defined by (5), according to the approach adopted here. There are some research work concerning these issues but good understanding of them has not yet been achieved.

Previous experimental and numerical results have shown that the form of the excitation profile has a profound effect on the topology of the error surface (Cherkaeva and Tripp, 1996; Figueroa and Seleghim, 2001). In the case of classical Dirac excitation, the application of a punctual voltage and grounding of the rest of the boundary is ineffective because most of the sensing energy tends to flow out of the measurement volume through electrodes near the excitation site. This implies in an intrinsic lack of sensitivity, especially in the central regions of the domain, not adequately illuminated by the sensing energy, and is directly related to the plateau surrounding the solution minimum. A possible solution to this consists in constructing an excitation profile that induces the electric sensing field to cross the measurement volume. This can be accomplished by applying a voltage profile that varies gradually from a maximum to a minimum value occurring at opposite points: a current line will probably not be deflected because the regions neighboring its injection site have similar potentials. Different excitation profiles can be generated according to this principle. Starting from the classic Dirac profile (Fig.3a) a pyramidal distribution can be generated by varying the potential from its maximum value of 10V, at the center of the upper boundary surface, to a minimum value of zero occurring on the median segment of the opposed side (Fig.3b). Since the pyramidal distribution presents some practical difficulties to be put in practice, an alternative profile which roughly follows the principle mentioned above can be generated by translating the excitation site along the upper median segment. This profile will be denoted ridge-Dirac and is shown in Fig.3c.

These excitation profiles were used in the numerical simulations to generate error surfaces e(m,n), showed in Figures 4 and 5, for approximated contrast defined according to the previous example, indicated in Fig.1. In order to provide a good comparison and to have a precise idea of the kind of improvement that can be achieved through these optimized excitation profiles, the former results (Fig.2) corresponding to classical Dirac excitation were also included in Figures 4 and 5. Considering the implementation of a numerical minimization procedure capable of autonomously locating the solution global minimum, that is, of refining the approximate permittivity distribution so that s_approx® s_actual, some main aspects concerning the nature of these error surfaces must be stressed out. First, the distinction between the solution minimum and the rest of the error surface has been significantly enhanced, what contributes to extending its attraction region. Second, the topology within these dominance regions have also been improved in the sense of being better distributed and assuring good convergence conditions for minimization methods that take the local inclination into account. It is also clear that some of the negative features concerning the convergence, such as the contour ridges, have been altered in a way that it is possible to find a minimization trajectory leading uniformly to the true minimum.

All these error surfaces can be combined into the global error mfunction defined in Eq.11 and, by doing so, some of the negative features of its topology are improved. In what concerns the success of a numerical minimization procedure, probably the most important improvement was obtained on the plateau surrounding the solution minimum, within which the convergence rate is significantly reduced. As it can be observed in Fig.6 this flat region has been averaged out and a uniform inclination around the solution was obtained. A much better enhancement of the inclination can be obtained through the error functional defined in (12), although multiple minima are introduced as it can be observed in Fig.7.

Genetic Algorithms and its Specialization to the EIT Problem

A Genetic Algorithm (GA) is a stochastic search algorithm based on the concepts of evolutionary theory (Goldberg, 1989; Holland, 1975). More specifically, solution candidates of the optimization problem (individuals) compete among themselves for the opportunity of transmitting their characteristics to new generations of solution candidates generated from the old ones (reproduction). The best fitted individuals have better chances of reproducing, i.e. those for which the optimization function evaluates to more optimized values are more likely to reproduce (selection), so that good characteristics tend to be preserved in new generations (elitism). Reproduction is also affected by random changes occurring at a certain probability (mutation), which enables the emergence of new characteristics, not initially present in the parent generations, and their exploration through the survival or decline of the mutant individuals. The application of these concepts into a particular optimization problem requires a special formalism which will be described in the sequel for the EIT problem treated in this work.

Considering the error function defined in (8) and the discretization proposed in (13) and (14), it is possible to establish that

The excitation profile being fixed and imposed externally, the error depends only on the contrast vector which must be refined to minimize (16). In other words, Eq.16 will be chosen as the fitness function and the corresponding individuals (or chromosomes) will be defined as

where N represents the total number of nodal points in W (N = N_x.N_y.N_z).The i-th contrast value in (17) corresponds to a gene of a particular individual and is related to each nodal contrast s_i,j,k though the following expression

in which the indexes i, j and k vary according to (15). Depending on the representation of the genes the GA is referred to be binary-coded, such as proposed by Holland (1975), or floating point-coded as employed by Michalewicz et al. (1994) and Cho et al. (2001) for dealing with EIT problems. Real-coding is more suited to multidimensional problems with non-trivial restriction and will be adopted in this work (Michalewicz, 1996).

A generation is defined as a collection of M individuals

which will reproduce to form the succeeding generation, according to the corresponding values of the fitness function so that the best fit ones have better chances of reproducing. The reproduction process is accomplished by combining the genes of two or more parent individuals. This can be done by merging randomly defined genes segments (geometric crossover: Radicliffe, 1990; Spears and De Jong, 1991) or by summing them up according to specific weights (arithmetic crossover: Davis, 1991; Gen and Cheng, 1997; Mühlenbein and Schlierkamp-Voosen, 1993).

In this work we adopted an extended arithmetic reproduction strategy with P parents combining to generate one descendent according to the following formula:

where the weights l_i vary randomly restricted to

and d is a positive parameter governing the extension outside the hyper-polygon formed by the parents in which descendents can be generated. This procedure has an important advantage over standard weighted averaging (linear crossover) defined by

because in this last the descendents are no homogeneously distributed over the parents hyper-polygon. To better illustrate this point consider the 2D example in which the parents, given by

₁ = (1.0,0.2), ₂ = (0.2,1.0) and ₃ = (1.3,1.2), combine to generate 50000 descendents according to (20) and (22), with random weights satisfying (21) and d=1.3. The following figure shows the polygon formed by these three parents (triangle vertices) and the corresponding descendents (dots).

The reproduction scheme produces a steady convergence to the optimal solution but suffers from some drawbacks as, for instance, the trapping by local minima or the dramatic convergence decrease due to plateaus in the optimization hyper-surface. To avoid this, from time to time, at a given probability, some of the reproduction rules are intentionally broken in a way that the generated descendents may acquire new features not necessarily present in the parent individuals. This process, called mutation, also assures that no point in the optimization hyper-surface has a null probability of being considered as a solution candidate, which implies that reaching the global optimal is at least possible, independently of initial guesses or how complex the pathology of the problem is. Mutation can be implemented according to different approaches depending on the coding of the chromosome (binary or real coded). The mutation strategy adopted in this work consists in selecting a random position in the chromosome and changing the corresponding gene by another number randomly chosen between predefined lower and upper limits (Cho et al, 2001). These limits correspond to s_min = 1 and s_max = 80 representing respectively continuous air phase and continuous distillated water phase.

A very important advantage of genetic algorithms in what regards their use in the solution of inverse ill-posed problems is the possibility of introducing additional restrictions representing information known a priori. As already mentioned, without this, the problem becomes extremely sensitive to experimental and round off errors, and probably unsolvable under realistic experimental conditions. In our case, it is possible to calculate from the signals delivered by a direct imaging probe, and appropriate calibration models, the instantaneous volumetric void fraction of the flow (a), in addition to the initial guess contrast distribution within the sensing region. In mathematical terms this can be expressed as

where P stands for the calibration model. For the sake of simplicity we adopted P(x,y,z) º 1 so that the void fraction becomes simply the volume integral of the internal contrast distribution. Thus, reproduction and mutation may be followed by an additional restriction operator which assures that all descendents comply with (23). This can be done by simply rejecting the unviable descendents (which is computationally inefficient) or by modification of the reproduction and mutation algorithms (which is not always feasible). We adopted a strategy based on reparation of the unviable descendents as described in Davis (1991), where an extensive survey on such procedures can also be found.

The last relevant genetic operation is the selection of the best fit that will reproduce to form the next generation. There are several selection strategies, each with a number of possible variations depending on specific requirements or computational performance criteria. The most common strategy, called stochastic, consists in assigning selection probabilities according to the fitness function of each individual, in a way that the best fit ones have greater chances of being selected for reproduction. In this work we adopted the so called deterministic selection, which consists in ranking the individuals based on their fitness value and selecting a predefined number of the best ones (elitism). Both parents and descendents were considered in the ranking process to enhance crossover possibilities and to assure that a succeeding generation is never worst than the previous one.

The basic implementation of our genetic algorithm can be summarized in the following steps:

1. An initial generation of chromosomes is constructed with random contrast values in (19), each complying with an overall structure given by the signals of the direct imaging probe (initial guess).

2. A priori information regarding the void fraction is enforced by reparation of the initial generation so that all individuals satisfy (23). These constitutes the first parent generation.

3. The fitness function of each chromosome is evaluated: Eq.1 is solved with the proper excitation condition given by (2), the corresponding response is calculated by means of (3) and its difference with respect to the correct measured response is assessed by (8).

4. The best fit chromosomes are selected for reproduction, i.e. only those with the smaller values of (16), according to the rule given by (20). A small percentage of descendents suffer mutation at random positions in the chromosome.

5. A priori information regarding the void fraction is enforced by reparation of the descendent chromosomes so that all individuals satisfy (23). The offspring generation is thus created.

6. The fitness function of each offspring chromosome is evaluated as described in step 3 above. Convergence is achieved if the minimum fitness value is smaller than a predefined target and the solution is the associated chromosome. The algorithm is equally stopped if the total number of generations exceeds an allowable maximum.

7. Both parent and offspring generations are ranked according to the corresponding fitness function and the best ones are selected for reproduction and mutation as describe in the step 5 above.

8. Iterate steps 5 to 7 until convergence or the maximum number of generations is achieved.

It is important to stress that the general performance of any genetic algorithm depends on a fine tuning of intrinsic parameters such as the number of individuals in the parent and offspring generations, number of combining parents and overextension of the corresponding hyper-polygon (respectively parameters P and d in Eq.20), mutation percentages and so forth. Intrinsic procedures such as the reproduction algorithm or the reparation strategy have also a strong influence on the overall performance and must be considered in detail. These points were cleared in this work by extensive numerical tests in which intrinsic parameters and procedures were set so guarantee the best convergence of the test case depicted in Fig.1. Whether this tuning is case-dependent or not and what would be a case-independent tuning implementation of a genetic algorithm is beyond the scope of this work and certainly represents an interesting and still open research topic.

Numerical Tests and Results

The error surface analysis described in previous section revealed important problematic characteristics of the optimization surfaces regarding the convergence of an optimization method to the global minimum. These characteristics, which define the so-called pathology of the optimization surface, were closely considered in the implementation and tuning of the genetic algorithm described in the previous section. Thus, two series of numerical tests were performed to demonstrate the ability of achieving convergence to the correct solution, that is, the contrast distribution shown in Fig.1. In the first series of tests only the position of the inclusion is unknown, being its shape, size and contrast values known by hypothesis. In the second series of tests all of this is unknown, except that there is an inclusion in a known search region contained in the sensing domain. The most important difference between these tests regards the dimensionality of the associated optimization problem: at maximum 3D in the first case while in the second case the dimension is equal to the number of discretization nodes defining the search region. In the first series of tests all the three excitation profiles given in Fig.3 were considered in order to demonstrate the influence of an improved excitation strategy. The corresponding results will be shown in the sequel.

Adopting the indices m, n and o to define the position of the frontal lower left corner of the prospective inclusion, as indicated in Fig.1, the error function (16) can be written as

Allowing the prospective inclusion to move along the transversal plane used to construct the optimization surfaces of Fig.4 (vertical search), the corresponding optimization problem can be mathematically defined as:

Note that in (25) the total number of possible positions that can be occupied by the inclusion is N_x.N_y = 400, which corresponds to the number of times that the fitness function must be evaluated in an exhaustive search for the minimum. In other words, the genetic algorithm or any other optimization method should be capable of finding the solution in much less than this in order to be advantageous. Similar tests were performed allowing the inclusion to move along the longitudinal plane (horizontal search, Fig.5), that is

in which case an exhaustive search would require N_x.N_z = 1200 evaluations of the fitness function, and also inside the whole sensing domain (3D search), that is

where the total number of evaluations of the fitness function is N_x.N_y.N_z = 24000.

The solution of these problems, i.e. (25), (26) and (27), was accomplished through the procedure outlined in the previous section, with the following specific parameters and procedures:

• Parent generation: 10 chromosomes randomly generated (the initial guess provided by the direct imaging probe is disregarded) and satisfying the internal contrast intrinsic limits (1 < g_i< 80)

• A priori restrictions: void fraction is intrinsically satisfied because the size of the prospective inclusion is the same of the real inclusion, so there is no need for reparation

• Reproducing population: the first and the second best fit chromosomes combining to generate 4 descendents.

• Mutation: two individuals in each generation are forced to mutate; only a single gene gi is replaced by a random value satisfying the internal contrast intrinsic limits

• Selection: 10 best fit chromosomes are selected out of 14 (4 descendents plus 10 parents) according to their fitness function values

• Target values: 10^-14, which represents to the precision of the machine, or 300 generations

The corresponding results are shown in Figures 9 to 11, being one for each type of search (vertical, horizontal or 3D). Each figure has three graphs for three different excitation profiles (classical Dirac, pyramidal and ridge Dirac). The fitness values of the best fit and of the worst fit individuals are plotted in each graph in order to assess the evolution of the population. These convergence results are summarized in table 1.

Several interesting features of the problem can be understood from these results. First of all, in all three cases, the pyramidal excitation profile produced the fastest convergence to the correct solution. The ridge Dirac excitation profile also gives reasonable results, with the advantage of being considerably simpler to implement in practice. The worst results were obtained from the classical Dirac excitation profile, having even failed to converge after 1210 evaluations of the fitness function in the 3D problem (Fig.11 on top). It is also interesting to notice that convergence is achieved in proportionally less iterations for the 3D search cases, i.e. those for which an exhaustive search is probably prohibitive in practical applications. This confirms the ability of genetic algorithms in properly dealing with high dimensionality optimization tasks and sustains that they constitutes a good approach for solving EIT problems.

Another important feature regards the so-called premature convergence, namely when the best and the worst fit chromosomes, and by consequence all the parent population, have virtually the same fitness value (the difference is smaller than the precision of the machine). This happens when all the reproducing chromosomes are located on the plateau around the global minimum of the optimization surfaces (see Figs. 4 and 5). Although an optimization method based exclusively on the local inclination would certainly fail, genetic algorithms rely on mutation to resume convergence after some generations, as can be observed for instance in Figs. 9, 10 and 11. How many is "some" generations depends strongly on the pathology of the optimization surface: while premature convergence was never observed in the simulations with pyramidal excitation it was observed in all tests with the classical Dirac excitation profile. Also, comparing the classical and the ridge Dirac convergence curves (both exhibit premature convergence) it is possible to observe that the latter in all tests needed fewer generations to resume convergence.

As mentioned above, a second series of tests was carried out to investigate the performance of the genetic algorithm in the solution of a high dimensionality optimization problem. Aiming at limiting this dimension to a reasonable value, from a computational point of view, the inclusion was retrieved from a three-dimensional search region contained in the sensing domain. This search region, shown in Figure 12, corresponds to a 3 ´ 3 ´ 3 cube (in mesh steps) and comprises 64 nodes whose corresponding contrast values must to be determined. In this case, not only the dimensionality of the optimization problem is considerably higher, but also all the 64 optimization variables may assume any values between s_min = 1 and s_max = 80. In mathematical terms this problem can be formulated according to

The inclusion of supplementary a priori information is of crucial importance for the problem to be solvable in a reasonable number of iterations. Thus, in addition to the void fraction as expressed by (23), another a priori parameter will be considered: the symmetry degree (c) of the contrast distribution within the search volume (as proposed by Cohen, 1995), which is given by

To enforce the symmetry degree it is necessary to adapt the genetic operators so that every descendent comply with (29), as was done for the void fraction in the examples above, or the fitness function can be modified to include a penalty term such as

which will penalize the individuals not complying with (29).

In this second series of tests, the following situations were simulated in order to emphasize the importance of introducing additional a priori information: 1) no a priori restriction, 2) void fraction only, 3) symmetry degree only and 4) both void fraction and symmetry degree are enforced. Hence, the solution of problem defined by (28) was achieved with specific parameters and procedures as follows:

• Parent generation: 200 chromosomes in which the 64 values in the 3D search region defined in

  Fig.12 are randomly generated satisfying the internal contrast intrinsic limits (1

  < g

  _i

  < 80)

• Void fraction: enforced through reparation of the chromosomes: a) if a

  _actual > a

  _approx then randomly chosen genes g

  _k in (28) are increased until (23) is satisfied and b) if a

  _actual > a

  _approx then randomly chosen genes g

  _k in (28) are randomly lowered until (23) is satisfied

• Symmetry degree: enforced through penalization of the fitness function according to (30), that is, the penalty is proportional to the difference between c

  _actual and c

  _approx

• Reproducing population: two chromosomes are randomly chosen from the first 20 best fit parents to generate 20 descendents by arithmetic crossover (Eq.20 with d = 10) and 20 descendents by geometric crossover

• Mutation: each mutant chromosome has 20% of its genes randomly replaced by s

  _min or s

  _max, with equal chances. Each generation has 60 mutants.

• Selection: the 200 parents and the 100 descendents are altogether ranked, the first 200 chromosomes are selected for the succeeding parent population

• Target value: 10

  ^-4 or 100 generations (approximately 10000 evaluations of the fitness function).

Pyramidal excitation was used in all cases aiming to enhance convergence. A Pentium III based computer (1.2GHz single processor with 1GB memory) was used for the calculations. In this machine, each evaluation of the fitness function takes approximately 4 seconds, mainly expended to solve the N_x.N_y.N_z = 24000 linear equations given by (13). A modified version of the conjugate gradient method given in Press et al. (1992) was used to solve these equations (the problem matrix was assembled directly on the storage vectors). Other significant time consuming tasks are the generation of random numbers and the multiplications implied in the extended arithmetic reproduction strategy (Eq.20).

Convergence results and the corresponding chromosomes are shown in Figure 13. It is clear from these graphs that the simultaneous enforcement of void fraction and symmetry restrictions significantly enhances the quality of the reconstructed chromosome (cf the correct one in Fig.12). After the same amount of evaluations of the fitness function the other simulations produced poor results, particularly when no a priori restriction is enforced. It is also interesting to notice that, in the void fraction only case, the solution has probably been trapped by a local minimum surrounded by a hyper-plateau such as the one in Fig.5 for classical Dirac excitation.

Conclusions

The application of electrical impedance tomography techniques in industrial processes is of great interest, particularly for those involving multiphase flow equipments due to the strong relation between their global performance and the characteristic organization of the flowing phases (flow regimes). Despite this, from an industrial point of view, there are still several issues to be considered, among which probably the most important one refers to consistency aspects of the final images.

The approach adopted in this work consists in formulating the reconstruction problem (as defined by Eq.5) in terms of an error functional expressing the difference between an approximated and the actual contrast distribution (image), from which one has access only indirectly through measurements of the boundary charge distribution that results from given boundary excitation conditions.

Numerical simulations have been carried out aiming to demonstrate the feasibility of our approach. The sensing volume was a 1:1:3 parallelepiped with voltage profiles imposed on its longitudinal boundaries (excitation) and no-flux condition imposed on the transversal boundaries. The flowing two-phase mixture corresponds to distillated water (s = 80) with a small three-dimensional air inclusion (s = 1) placed as specified in Fig.1, which is obviously not a realistic situation (a cubic bubble) but represents a much more difficult measurement condition and, consequently, a more severe test to the proposed formulation of the reconstruction problem. Approximated contrast distributions are generated by translations of a prospective inclusion and the resulting electrical charge distributions, calculated with (3) and the appropriate boundary condition, are used in (8) to construct error surfaces (Figs. 4, 5 and 6).

It is clear from these results that a close match between the actual and the approximated contrast distributions is associated with a pronounced minimum on the error surface, which indicates that, positively, a formulation as proposed in (8) allows the solution of the inverse problem (5) and to reconstruct the internal constitution of the sensed medium. This approach has several important advantages: the formulation is so that a true 3-dimensional reconstruction is achieved with no restrictive additional assumptions and free of averaging effects for sufficiently extended and dense mashing of the measurement volume. Another important advantage concerns the straightforwardness of introducing additional useful information for the reconstruction procedure (regularization, redundancy, etc.). A good example of the importance of this last feature is the development of multimodal tomographic techniques, in which objects faintly observable with one field presumably could be observed with other type of sensing field.

Some drawbacks are also evident and will certainly demand research efforts to overcome them. Probably the most important one is that the procedure proposed here is intrinsically iterative, thus, in principle, not well suited for real time implementation. In spite of the difficulties that may arise from this, the fact that direct imaging techniques may provide very good initial approximations of the actual image (thus demanding a few refinement iterations) and the exponentially growing computational power of most common platforms (processing speed, parallelization, etc.) will, in our belief, make this an irrelevant issue.

Acknowledgements

Support to this work was provided by FAPESP through grants 98/12921-1 and 99/02821-2.

Paper accepted October, 2005.

Technical Editor: Aristeu da Silveira Neto.

• Allers A.; F. Santosa, 1991, "Stability and resolution analysis of a linearized problem in electrical impedance tomography", Inverse Problems, vol. 7, pp. 51533.
• Borcea L., 2003, "Electrical impedance tomography", Inverse Problems, vol. 18, no. 6, pp. R99-R136.
• Calderon A.P., 1980, "On an inverse boundary value problem", Seminar on Numerical Analysis and its applications to Continuum Physics Brazilian Society of Mathematics, pp. 6573.
• Cherkaeva E. and A. Tripp, 1996, "Inverse conductivity problem for noisy measurements", Inverse Problems, vol. 12, pp. 869883.
• Cho K.H., S. Kim and Y.J. Lee, 2001, "Impedance imaging of two-phase flow field with mesh grouping method", Nuclear Engineering and Design, vol. 204(1-3), pp. 57- 67.
• Cohen L., 1995, "Time Frequency Analysis", Prentice-Hall, Englewood Cliffs, 320p.
• Davis L., 1991, "Handbook of Genetic Algorithms", Van Nostrand Reinhold, 385p.
• Dobson D.C., 1992, Estimates on resolution and stabilization for the linearized inverse conductivity problem, Inverse Problems, vol. 8, pp. 8 71881.
• Engl H.W., M. Hanke and A. Neubauer, 1996, "Regularization of Inverse Problems (Mathematics and its Applications)", Kluwer Academic Publishers, vol. 375, 321p.
• Figueroa T.P. and P. Seleghim Jr., 2001, "Sensitivity Analysis of Different Sensing Strategies for Electrical Impedance Imaging of Two-Phase Flows", Journal of Electronic Imaging, vol. 10, no. 3, pp. 641-645.
• Gen M. and R. Cheng, 1997, "Genetic Algorithms & Engineering Design", Wiley Interscience Publication, 432p.
• Goldberg D.E., 1989, "Genetic algorithms in search, optimization, and machine learning", New York: Addison-Wesley, 432p.
• Hadamard J., 1902, "Sur les problèmes aux dérivées partielles et leur signification physique", Bulletins of the University of Princeton, vol. 13, pp.49-52.
• Holland J.H., 1975, "Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology", Control, and Artificial Intelligence, The University of Michigan Press, 211p.
• Hsiao C.T., G. Chahine and N. Gumerov, 2001. "Application of a Hybrid Genetic/Powell Algorithm and a Boundary Element Method to Electrical Impedance Tomography", Journal of Computational Physics, vol. 173, pp.433454.
• Knowles I., 1998, "A variational algorithm for electrical impedance tomography", Inverse Problems, vol.14, pp. 15131525.
• Kohn R.V. and M. Vogelius, 1984, "Determining conductivity by boundary measurements", Communications on Pure and Applied Mathematics, vol. 37, pp. 113123.
• Kohn R.V. and M. Vogelius, 1985, "Determining conductivity by boundary measurements: 2. interior results", Communications on Pure and Applied Mathematics, vol. 38, pp. 643867.
• Lányi S., 1998, "Analysis of linearity errors of inverse capacitance position sensors", Measurement Science and Technology, Vol. 9, pp.1757-1764.
• Michalewicz Z., 1996, "Genetic Algorithms + Data Structures = Evolution Programs", New York: Springer-Verlag Berlin Heidelberg, 387p.
• Michalewicz Z., T. Logan and S. Swaminathan, 1994, "Evolutionary operators for continuous convex parameter spaces", Proceedings of the 3rd Annual Conference on Evolutionary Programming, A.V. Sebald and L.J. Fogel (editors), World Scientific Publishing, pp.84-97.
• Mühlenbein H. and D. Schlierkamp-Voosen, 1993, "Predictive Models for the breeder genetic algorithm I. continuous parameter optimization", Evolutionary Computation, vol. 1, pp. 25-49.
• Nachman A.I., 1988, "Reconstructions from boundary measurements", Annals of Mathematics, vol. 128, pp. 531576.
• Nachman A.I., 1996, "Global uniqueness for a two-dimensional inverse boundary problem", Annals of Mathematics, vol. 143, pp. 7196.
• Press W.H., B.P. Flannery., S.A. Teukolsky W.T. and Vetterling W. T., 1992, "Numerical Recipes in Fortran", Cambridge University Press, 966p.
• Radicliffe N., 1990, "Genetic neural networks on MIMD Computers", Ph.D. thesis, University of Edinburgh, UK, 90p.
• Rolnik V.P. and P. Seleghim Jr., 2002, "On site calibration of a phase fraction meter by an inverse technique", Journal of the Brazilian Society of Mechanical Sciences and Engineering, v.24, n.4, p.1 14.
• Seleghim Jr. P. and E. Hervieu, 1998, "Direct imaging of horizontal gas-liquid flows". Measurement Science & Technology, vol.9, no.8, pp.1492-1500.
• Seleghim Jr. P. and F. Schiavon, 2001, "On-site calibration of a phase fraction meter by an inverse technique", Proceedings of the 4th International Conference on Multiphase Flow, New Orleans USA (2001).
• Seleghim Jr. P. and F.E. Milioli, 2001, "Improving the determination of bubble size histograms by employing wavelet de-noising techniques", Powder Technology, vol. 115, pp.114-123.
• Spears W. and K. De Jong, 1991, "On the virtues of parameterized uniform crossover", Proceedings of the Fourth International Conference on Genetic Algorithms, San Mateo USA, pp. 230 236.

Publication Dates

• Publication in this collection
  08 Oct 2007
• Date of issue
  Dec 2006

History

• Accepted
  Oct 2005
• Received
  Oct 2005