Zbigniew Michalewicz

Marc Schoenauer

Evolutionary Algorithms for Constrained Parameter Optimization Problems

Evolutionary computation techniques have received a lot of attention regarding their potential as optimization techniques for complex numerical functions. However, they have not produced a signicant breakthrough in the area of nonlinear programming due to the fact that they have not addressed the issue of constraints in a systematic way. Only recently several methods have been proposed for handling nonlinear constraints by evolutionary algorithms for numerical optimization problems; however, these methods have several drawbacks and the experimental results on many test cases have been disappointing.

In this paper we (1) discuss diculties connected with solving the general nonlinear programming problem, (2) survey several approaches which have emerged in the evolutionary computation community, and (3) provide a set of eleven interesting test cases, which may serve as a handy reference for future methods.

Introduction

The general nonlinear programming (NLP) problem is to nd x so as to optimize f(x), x = (x 1 ; : : : ; x n ) 2 I R n , where x 2 F S. The objective function f is dened on the search space S I R n and the set F S denes the feasible region. Usually, the search space S is dened as a n-dimensional rectangle in I R n (domains of variables dened by their lower and upper bounds): l(i) x i u(i); 1 i n, whereas the feasible region F S is dened by a set of m additional constraints (m 0): g j (x) 0, for j = 1; : : : ; q, and h j (x) = 0, for j = q + 1; : : : ; m. At any point x 2 F, the constraints g k that satisfy g k (x) = 0 are called the active constraints at x. By extension, equality constraints h j are also called active at all points of S.

The NLP problem, in general, is intractable. If the objective function f and functions g j 's and h j 's expressing constraints are arbitrary, then there is little choice apart from methods based on exhaustive search. This is the reason for identifying several special cases of the NLP. One of the studied cases is the one where all functions g j 's and h j 's are linear; such problems are called linearly constrained optimization problems. If, additionally, the objective function f is polynomial at most quadratic, the problem is called quadratic programming problem. The best known special case of quadratic programming is where the objective function f is linear as well; this problem is called linear programming problem. There is also an important special case, called unconstrained optimization, where there are no constraints at all, i.e., m = 0 and F = S.

One of the main diculties in solving NLP is the problem of local optima; a feasible point x0 2 F is a local optimum for the objective function f i there exist > 0 such that for all x in the -neighborhood of x0 in F, f(x) > f( x0 ). 1 Local optima just satisfy the mathematical requirements on the derivatives of the functions and many optimization techniques based on gradient methods aim at local optimization only.

In general, it is impossible to develop a deterministic method for the NLP in the global optimization category, which would be better than the exhaustive search. As stated by [START_REF] Gregory | Nonlinear Programming FAQ[END_REF] \It's unrealistic to expect to nd one general NLP code that's going to work for every kind of nonlinear model. Instead, you should try to select a code that ts the problem you are solving. If your problem doesn't t in any category except `general', or if you insist on a globally optimal solution (except when there is no chance of encountering multiple local optima), you should be prepared to have to use a method that boils down to exhaustive search, i.e., you have an intractable problem."

Evolutionary algorithms are global methods, which aim at complex objective functions (e.g., non dierentiable or discontinuous). However, most research on applications of evolutionary computation techniques to nonlinear programming problems has been concerned with complex objective functions but no constraints (F = S). Many of the test functions used by various researchers during the last 20 years did not include any constraints (apart from specied domains 1 For minimization problems. 2 of variables); this was the case with ve test functions F1{F5 proposed by De Jong (1975), as well as with many other test cases proposed since then [START_REF] Eshelman | Real-coded genetic algorithms and interval-schemata[END_REF][START_REF] Fogel | On the eectiveness of crossover in simulated evolutionary optimization[END_REF][START_REF] Wright | Genetic algorithms for real parameter optimization[END_REF]. Only recently several approaches have emerged which aim at general nonlinear programming problems.

It seems worthwhile to extend evolutionary techniques by some constraint-handling methods. It is also important to gain some experimental feedback on merits and drawbacks of these methods. The ultimate goal is to create a useful evolutionary toolbox for the general NLP problem!

The paper is organized as follows. The following section surveys briey several operators which has been developed for evolutionary algorithms with oating-point representation. Section 3 presents several constraint-handling techniques for numerical optimization problems which have emerged in evolutionary computation techniques over the years. Most methods are enlightened by an example of how they behave on a specic NLP problem. Section 4 summarizes all methods and indicates some directions for future research.

2 Numerical Optimization and Operators Most evolutionary algorithms for numerical optimization problems use vectors of oating point numbers for their chromosomal representations. For such representations, many operators have been proposed during the last 30 years. Several constraint-handling methods are discussed in the following section; some of these methods (e.g., Genocop, section 3.1.1) require a set of specic operators. For that reason, we discuss briey these operators in turn.

Mutation operators

The most popular mutation operator is Gaussian mutation, which modied all components of the solution vector x = hx 1 ; : : :; x n i by adding a random noise: xt+1 = xt + N(0;), where N(0;) is a vector of independent random Gaussian numbers with a mean of zero and standard deviations . Such a mutation is used in evolution strategies (B ack et al., 1991) and evolutionary programming [START_REF] Fogel | Evolutionary Computation. Toward a New Philosophy of Machine Intelligence[END_REF]. (One of the historical dierences between these techniques lies in adjusting vector of standard deviations ).

Other types of mutations include non-uniform mutation, where

x t+1 k = ( x t k + 4(t; r(k) x k ) if a random binary digit is 0 x t k 4(t; x k l(k)) if a random binary digit is 1 for k = 1; : : : ; n. The function 4(t; y) returns a value in the range [0; y] such that the probability of 4(t; y) being close to 0 increases as t increases (t is the generation number). This property causes this operator to search the space uniformly initially (when t is small), and very locally at later stages. In experiments reported by Michalewicz et al. (1994), the following function was used: 4(t; y) = y r (1 t T ) b ; where r is a random number from [0::1], T is the maximal generation number, and b is a system parameter determining the degree of non{uniformity.

It is also possible to experiment with a uniform mutation, which changes a single component of the solution vector; e.g., if xt = (x 1 ; : : :; x k ; : : : ; x n ), then xt+1 = (x 1 ; : : : ; x 0 k ; : : :; x n ), where x 0 k is a random value (uniform probability distribution) from the domain of variable x k . A special case of uniform mutation is boundary mutation, where x 0 k is either left of right boundary of the domain of x k .

Crossover operators

There are several interesting types of crossover operators. The rst family is directly designed from the well-known bitstring crossover operators (known also as \discrete recombination" operators): 1-point, 2-point and uniform crossover, except that the \genes" exchanged between the two parents are real-valued coordinates. For instance, uniform crossover, works as follows.

Two parents, x = hx 1 ; : : :; x n i and ỹ = hy 1 ; : : :; y n i, produce an ospring, z = hz 1 ; : : :; z n i, where z i = x i or z i = y i with equal probability for all i = 1; : : : ; n (the second ospring is created by reversing decisions for all components). Some multi-parents versions of these operators have also been designed. M uhlenbein and [START_REF] Ack | Gene poool recombination for the breeder genetic algorithm[END_REF] investigated the properties of a recombination operator called gene pool recombination (default recombination mechanism in evolution strategies [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF], where the genes are randomly picked from the gene pool dened by the selected parents. A similar multiparent crossover was also investigated also by [START_REF] Eiben | Genetic algorithms with multi-parent recombination[END_REF] in the context of combinatorial optimization.

Another possibility is arithmetical crossover. Here two vectors, x and ỹ produce two ospring, p and q, which are linear combinations of their parents, i.e., p = a x + (1 a) ỹ and q = (1 a) x + a ỹ. This operator has also been called guaranteed average crossover [START_REF] Davis | Adapting operator probabilities in genetic algorithms[END_REF] (when a = 1=2), and intermediate crossover [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF].

It is possible to generalize arithmetical crossover into multi-parent operator; parents x1 ; x2 ; : : : ; xr may produce a family of ospring: ỹ = a 1 x1 + a 2 x2 + : : : + a r xr , for any choice of (a i ) 2 [0; 1] r such that a 1 + a 2 + : : : + a r = 1.

The simplex crossover operator proposed by [START_REF] Renders | Hybridizing genetic algorithms with hill-climbing methods for global optimization: Two possible ways[END_REF] for numerical optimization problems can also be viewed as another generalization of arithmetical crossover; this crossover operator involves computing the centroid of group of parents and moving from the worst individual beyond the centroid point. More precisely, the operator selects k > 2 parents (set J), determines the best and the worst individual within the selected group ( b and w, respectively), computes the centroid c of the selected group with removed worst individual: c = P xi 2J w xi =(k 1), and computes the `reected point' ỹ (i.e., the ospring) obtained from the worst one: ỹ = c + (c w). (This operator includes also a few extra options, which are based on the outcome of comparisons between the values of f(ỹ), f(w), and f( b); for details, see [START_REF] Renders | Hybridizing genetic algorithms with hill-climbing methods for global optimization: Two possible ways[END_REF].

Some of these ideas were embodied earlier in the evolutionary procedure called scatter search [START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF]. The process generates initial populations by screening good solutions produced by heuristics. The points used as parents are then joined by linear combinations with contextdependent weights, where such combinations may apply to multiple parents simultaneously. The linear combination operators are further modied by adaptive rounding processes to handle components required to take discrete values. (The vectors operated on may contain both real and integer components, as opposed to strictly binary components.) Finally, preferred outcomes are selected and again subjected to heuristics, whereupon the process repeats. The approach has been found useful for mixed integer and combinatorial optimization problems. (For background and recent developments see, e.g., [START_REF] Glover | Genetic algorithms and scatter search: Unsuspected potentials[END_REF]).

An interesting variation along this line is the heuristic crossover operator proposed by [START_REF] Wright | Genetic algorithms for real parameter optimization[END_REF]; this crossover uses values of the objective function in determining the direction of the search, and it produces only one ospring. The operator generates a single ospring z from two parents, x and ỹ according to the following rule:

z = r (x ỹ) + x,
where r is a random number between 0 and 1, and the parent x is not worse than ỹ, i.e., f(x) f(ỹ) for maximization problems and f(x) f(ỹ) for minimization problems.

Recently, some results of experiments with a geometrical crossover operator were reported [START_REF] Michalewicz | A note on usefulness of geometrical crossover for numerical optimization problems[END_REF]. Geometrical crossover can only be applied to problems where each variable takes nonnegative values only, i.e., 0 l(i) x i u(i), so its use is quite restricted in comparison with other types of crossovers (e.g., arithmetical crossover). However, for many engineering problems, all problem variables are positive; moreover, it is always possible to replace variable x i 2 hl(i); u(i)i which can take negative values (i.e., where l(i) < 0) with a new variable y i = x i l(i).

The geometrical crossover operator takes two parents and produces a single ospring; for parents x and ỹ the ospring is z = h p x 1 y 1 ; : : : ; p x n y n i. Also it is possible to generalize this operator to include several parents: ỹ = h(x 1;1 ) 1 (x 1;2 ) 2 : : : (x 1;r ) r ; : : : ; (x n;1 ) 1 (x n;2 ) 2 : : : (x n;r ) r i, where 1 + : : : + r = 1 and 0 i 1 for all 1 i r.

Constraint-handling methods

During the last few years several methods were proposed for handling constraints by genetic algorithms for parameter optimization problems; this section discusses them in turn. It is useful to group these methods into four categories:

1. methods based on preserving feasibility of solutions, 2. methods based on penalty functions, 3. methods which make a clear distinction between feasible and infeasible solutions, and 4. other hybrid methods.

The subsequent four subsections 3.1{3.4 describe methods of these categories, respectively. Also, for most methods we provide with a test case and the results of the method on this test case. We do not dierentiate between minimization and maximization problems (we present test cases of both categories); note that min ff(x)g is equivalent to max f f(x)g, for x 2 F.

Methods based on preserving feasibility of solutions

There are two methods which fall into this category; we discuss them in turn.

Use of specialized operators

The Genocop (for GEnetic algorithm for Numerical Optimization of COnstrained Problems) system was developed by [START_REF] Michalewicz | Handling constraints in genetic algorithms[END_REF]. The idea behind the system is based on specialized operators which transform feasible individuals into feasible individuals, i.e., operators, which are closed on the feasible part F of the search space. The method assumes linear constraints only and a feasible starting point (or feasible initial population). Linear equations are used to eliminate some variables; they are replaced as a linear combination of remaining variables. Linear inequalities are updated accordingly. A closed set of operators maintains feasibility of solutions. For example, when a particular component x i of a solution vector x is mutated, the system determines its current domain dom(x i ) (which is a function of linear constraints and remaining values of the solution vector x) and the new value of x i is taken from this domain (either with at probability distribution for uniform mutation, or other probability distributions for non-uniform and boundary mutations). In any case the ospring solution vector is always feasible. Similarly, arithmetic crossover, ax + (1 a)ỹ, of two feasible solution vectors x and ỹ yields always a feasible solution (for 0 a 1) in convex search spaces (the system assumes linear constraints only which imply convexity of the feasible search space F). The heuristic crossover described in section 2.2 generates a single ospring z (note that an additional feasibility check is required in such a case) from two parents.

Genocop proved its usefulness for linearly constrained optimization problems; it gave surprisingly good performance on many test functions (Michalewicz et al., 1994;[START_REF] Michalewicz | Genetic Algorithms+Data Structures=Evolution Programs[END_REF]. For example, let us consider the following problem [START_REF] Floudas | A Collection of Test Problems for Constrained Global Optimization Algorithms[END_REF] to minimize a function: G1(x) = 5x 1 + 5x 2 + 5x 3 + 5x 4 5 P 4 i=1 x 2 i P 13 i=5 x i , subject to the following constraints:

2x 1 + 2x 2 + x 10 + x 11 10, 2x 1 + 2x 3 + x 10 + x 12 10, 2x 2 + 2x 3 + x 11 + x 12 10, 8x 1 + x 10 0, 8x 2 + x 11 0, 8x 3 + x 12 0, 2x 4 x 5 + x 10 0, 2x 6 x 7 + x 11 0, 2x 8 x 9 + x 12 0, and bounds 0 x i 1, i = 1; : : : ; 9, 0 x i 100, i = 10; 11; 12, 0 x 13 1.

The problem has 13 variables and 9 linear constraints; the function G1 is quadratic with its global minimum at x = (1; 1; 1; 1; 1; 1; 1; 1; 1; 3; 3; 3; 1), where G1(x ) = 15. Six (out of nine) constraints are active at the global optimum (all except the following three: 8x 1 + x 10 0, 8x 2 + x 11 0, 8x 3 + x 12 0). Genocop requires less than 1,000 generations to arrive at the global solution. Note that the linear constraints imply convexity of the feasible search space, which usually require smaller eort of the search process than non convex ones, where the feasible parts of the search space can be disjoint and irregular. Note also, that the method can be generalized to handle nonlinear constraints provided that the resulting feasible search space F is convex.

But the weakness of the method lies in its inability to deal with non convex search spaces (i.e, to deal with nonlinear constraints in general).

Searching the boundary of feasible region

One of the main reasons that lie behind diculties in locating the global solution is the inability of evolutionary systems to search precisely the boundary area between feasible and infeasible regions of the search space. This ability can be quite important in the case of optimization problems with nonlinear equality constraints or with active nonlinear constraints at the target optimum.

Some other heuristic methods recognized the need for searching areas close to the boundary of the feasible region. For example, one of the most recently developed approach for constrained optimization is strategic oscillation. Strategic oscillation was originally proposed in accompaniment with the strategy of scatter search [START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF], and more recently has been applied to a variety of problem settings in combinatorial and nonlinear optimization (see, for example, the review of [START_REF] Glover | Critical event tabu search for multidimensional knapsack problems[END_REF]). The approach is based on identifying a critical level, which represents a boundary between feasibility and infeasibility. The basic strategy is to approach and cross the feasibility boundary, by a design that is implemented either by adaptive penalties and inducements (which are progressively relaxed or tightened according to whether the current direction of search is to move deeper into a particular region or to move back toward the boundary) or by simply employing modied gradients or sub-gradients to progress in the desired direction.

Evolutionary computation techniques have a huge potential in incorporating specialized operators which search the boundary of feasible and infeasible regions in an ecient way. This might be quite signicant: it is a common situation for many constrained optimization problems that some constraints are active at the target global optimum. This optimum thus lies on the boundary of the feasible space. One the other hand, it is commonly acknowledged that restricting the size of the search space in evolutionary algorithms (as in most search algorithms) is generally benecial. Hence, it seems natural in the context of constrained optimization to restrict the search of the solution to the boundary of the feasible part of the space. We provide two examples of a such approach; for more details the reader is referred to [START_REF] Schoenauer | Evolutionary computation at the edge of feasibility[END_REF]).

An interesting constrained numerical optimization test case emerged recently; the problem [START_REF] Keane | Genetic Algorithms Digest[END_REF] is to maximize a function:

G2(x) = j P n i=1 cos 4 (x i ) 2 Q n i=1 cos 2 (x i ) pP n i=1 ix 2 i j;
subject to Q n i=1 x i 0:75 , P n i=1 x i 7:5n , and bounds 0 x i 10 for 1 i n.

Function G2 is nonlinear and its global maximum is unknown, lying somewhere near the origin. The problem has one nonlinear constraint and one linear constraint; the latter one is inactive around the origin and will be forgotten in the following. Some potential diculties of solving this test case are illustrated in gure 1, where infeasible points were assigned a value of zero. The interesting boundary between feasible and infeasible regions is dened by the equation x i = 0:75. It is a dicult problem, on which no standard method (be it deterministic or evolutionary) gave satisfactory results. As [START_REF] Keane | Genetic Algorithms Digest[END_REF] noted:

\I am currently using a parallel GA with 12 bit binary encoding, crossover, inversion, mutation, niche forming and a modied Fiacco-McCormick constraint penalty function to tackle this. For n = 20 I get values like 0.76 after 20,000 evaluations."

This test case was the rst one on which the idea of searching only the boundary was used [START_REF] Michalewicz | A note on usefulness of geometrical crossover for numerical optimization problems[END_REF]: due to the simple analytical formulation of the constraint, ad hoc specic initialization procedure and operators could be designed:

Initialization: Randomly choose a positive variable for x i , and use its inverse as a variable for x i+1 . The last variable is either 0.75 (when n is odd), or is multiplied by 0.75 (if n is even), so that the point lies on the boundary surface.

Crossover: The geometrical crossover is dened by (x i )(y i ) ! (x i y 1 i ), with randomly chosen in [0; 1] Figure 2 illustrates the possible ospring from two parents for all values of .

Mutation: Pick two variables randomly, multiply one by a random factor q and the other by 1 q (restrict q to respect the bounds on the variables). The geometrical crossover on the hyperboloid around its center: when parameter goes from 0 to 1, the ospring of two parents denes the line joining them on the plot.

The simple evolutionary algorithm described above gave outstanding results. For the case n = 20 the system reached the value of 0.80 in less than 4,000 generations (with population size of 30, probability of crossover p c = 1:0, and probability of mutation p m = 0:06) in all runs. The best value found (namely 0.803553) was better than the best values of any method discussed earlier, whereas the worst value found was 0.802964. Similarly, for n = 50, all results (in 30,000 generations) were better than 0.83 (with the best of 0.8331937).

It was interesting to note the importance of geometrical crossover. With xed population size (kept constant at 30), the higher values of probability of crossover p c , the better results of the system were observed. Similarly, the best mutation rates were relatively low (p m 0:06).

Since evolutionary methods for searching boundaries of feasible region constitute a new development in constrained optimization, let us illustrate this approach on another test case. The test problem was constructed for the case of sphere [START_REF] Michalewicz | A note on usefulness of geometrical crossover for numerical optimization problems[END_REF][START_REF] Schoenauer | Evolutionary computation at the edge of feasibility[END_REF]; the task is to maximize

G3(x) = ( p n) n Q n i=1 x i ,
where P n i=1 x 2 i = 1 and 0 x i 1 for 1 i n.

The function G3 has a global solution at (x 1 ; : : :; x n ) = ( 1p n ; : : :; 1 p n ) and the value of the function in this point is 1. The evolutionary algorithm uses the following components:

Initialization: Randomly generate n variables y i , calculate s = P n i=1 y 2 i , and initialize an individual (x i ) by x = y i =s for i 2 [1; n]. Crossover: The sphere crossover produces one ospring (z i ) from two parents (x i ) and (y i ) by: z i = q x 2 i + (1 )y 2 i i 2 [1; n], with randomly chosen in [0; 1]: Mutation: Similarly, the problem-specic mutation transforms (x i ) by selecting two indices i 6 = j and a random number p in h0; 1i, and setting:

x i ! p x i and x j ! q x j ; where q = s ( x i x j ) 2 (1 p 2 ) + 1:

The simple evolutionary algorithm described above gave very good results. For the case n = 20 the system reached the value of 0.99 in less than 6,000 generations (with population size of 30, probability of crossover p c = 1:0, and probability of mutation p m = 0:06) in all runs. The best value found in 10,000 generations was 0.999866.

Methods based on penalty functions

In Evolutionary Computation as in many other elds of optimization, most of the constrainthandling methods are based on the concept of (exterior) penalty functions, which penalize infeasible solutions, i.e. try to solve an unconstrained problem (on F) using the modied tness function: eval(x) = ( f(x); if x 2 F f(x) + penalty(x); otherwise; where penalty(x) is zero, if no violation occurs, and is positive, otherwise. Usually, the penalty function is based on the distance of a solution from the feasible region F, or on the eort to \repair" the solution, i.e., to force it into F. The former case is the most popular one; in many methods a set of functions f j (1 j m) is used to construct the penalty, where the function f j measures the violation of the j-th constraint in the following way: f j (x) = ( maxf0; g j (x)g; if 1 j q jh j (x)j; if q + 1 j m:

However, these methods dier in many important details, how the penalty function is designed and applied to infeasible solutions. In the following subsections we discuss all known methods in turn.

Method of static penalties

The method was proposed by [START_REF] Homaifar | Constrained optimization via genetic algorithms[END_REF]; it assumes that for every constraint we establish a family of intervals which determine appropriate penalty coecient. It works as follows:

for each constraint, create several (`) levels of violation, for each level of violation and for each constraint, create a penalty coecient R ij (i = 1; 2; : : : ; `, j = 1; 2; : : : ; m); higher levels of violation require larger values of this coecient.

start with a random population of individuals (feasible or infeasible), evolve the population; evaluate individuals using the following formula eval(x) = f(x) + P m j=1 R ij f 2 j (x). The weakness of the method is in the number of parameters. For m constraints the method requires m(2`+ 1) parameters in total: m parameters to establish number of intervals for each constraint, `parameters for each constraint, dening the boundaries of the intervals (levels of violation), and `parameters for each constraint representing the penalty coecients R ij . In particular, for m = 5 constraints and `= 4 levels of violation, we need to set 45 parameters! Clearly, the results are parameter dependent. It is quite likely that for a given problem there exist one optimal set of parameters for which the system returns feasible near-optimum solution, however, it might be quite hard to nd it.

A limited set of experiments reported by Michalewicz (1995a) indicates that the method can provide good results if violation levels and penalty coecients R ij are tuned to the problem. For example, [START_REF] Homaifar | Constrained optimization via genetic algorithms[END_REF] experimented with the following problem [START_REF] Himmelblau | Applied Nonlinear Programming[END_REF]. Minimize a function of 5 variables: G4(x) = 5:3578547x 2 3 + 0:8356891x 1 x 5 + 37:293239x 1 40792:141, subject to three double inequalities: 0 85:334407 + 0:0056858x 2 x 5 + 0:00026x 1 x 4 0:0022053x 3 x 5 92 90 80:51249 + 0:0071317x 2 x 5 + 0:0029955x 1 x 2 + 0:0021813x 2 3 110 20 9:300961 + 0:0047026x 3 x 5 + 0:0012547x 1 x 3 + 0:0019085x 3 x 4 25, and bounds: 78 x 1 102, 33 x 2 45, 27 x i 45 for i = 3; 4; 5.

The best solution obtained in 10 runs [START_REF] Homaifar | Constrained optimization via genetic algorithms[END_REF]) was x = (80:49; 35:07; 32:05; 40:33; 33:34) with G4(x) = 30005:7, whereas the optimum solution [START_REF] Himmelblau | Applied Nonlinear Programming[END_REF]) is x = (78:0; 33:0; 29:995; 45:0; 36:776), with G4(x ) = 30665:5. Two constraints (upper bound of the rst inequality and the lower bound of the third inequality) are active at the optimum.

However, the importance of proper values for penalty coecients should be emphasized. For example, [START_REF] Homaifar | Constrained optimization via genetic algorithms[END_REF] reported a value 50 for such penalty coecient for the violation of the lower boundary of the rst constraint and the value of 2.5 for the violation of the upper boundary of the rst constraints, without any comments on the source for these values.

Method of dynamic penalties

The method was proposed by [START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with gas[END_REF]. As opposed to the previous method, the authors assumed dynamic penalties. Individuals are evaluated (at the iteration t) by the following formula:

eval(x) = f(x) + (C t) P m j=1 f j (x),
where C, and are constants. A reasonable choice for these parameters, reported by [START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with gas[END_REF], is C = 0:5, = = 2. The method requires much smaller number (independent of the number of constraints) of parameters than the rst method. Also, instead of dening several levels of violation, the pressure on infeasible solutions is increased due to the (C t) component of the penalty term: towards the end of the process (for high values of the generation number t), this component assumes large values.

One of the problems [START_REF] Hock | Test Examples for Nonlinear Programming Codes[END_REF] experimented by [START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with gas[END_REF] was to minimize a function G5(x) = 3x 1 + 0:000001x 3 1 + 2x 2 + 0:000002=3x 3 2 subject to

x 4 x 3 + 0:55 0, x 3 x 4 + 0:55 0, 1000 sin( x 3 0:25) + 1000 sin( x 4 0:25) + 894:8 x 1 = 0 1000 sin(x 3 0:25) + 1000 sin(x 3 x 4 0:25) + 894:8 x 2 = 0 1000 sin(x 4 0:25) + 1000 sin(x 4 x 3 0:25) + 1294:8 = 0 0 x i 1200, i = 1; 2, and 0:55 x i 0:55, i = 3; 4.

The best known solution is

x = (679:9453; 1026:067; 0:1188764; 0:3962336), and G5(x ) = 5126:4981. The best result of experiments reported by [START_REF] Joines | On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with gas[END_REF] (many runs, various values of and ) gave a value of 5126.6653 of the objective function. No solution was fully feasible due to the three equality constraints, but the sum of the violated constraints was quite small (10 4 ). On the other hand, the method seems to provide too strong penalties: often the factor (C t) grows too fast to be useful. The system has little chances to escape from local optima: in most experiments reported by Michalewicz (1995a) the best individual was found in early generations. The method gave very good results for test cases where the objective functions were quadratic.

Method of annealing penalties

This method, called Genocop II, is also based on dynamic penalties and was described by [START_REF] Michalewicz | Evolutionary optimization of constrained problems[END_REF] and [START_REF] Michalewicz | Genetic Algorithms+Data Structures=Evolution Programs[END_REF]. Its modied version works as follows: divide all constraints into four subsets: linear equations, linear inequalities, nonlinear equations, and nonlinear inequalities, select a random single point as a starting point (the initial population consists of copies of this single individual). This initial point satises linear constraints, set the initial temperature = 0 , evolve the population using the following formula: eval(x; ) = f(x) + 1 2 P m j=1 f 2 j (x), if < f , stop, otherwise { decrease temperature , { the best solution serves as a starting point of the next iteration, { repeat the previous step of the algorithm. This is the only method we are aware of (except Genocop III, section 3.3.3) which distinguishes between linear and nonlinear constraints. The algorithm maintains feasibility of all linear constraints using a set of closed operators, which convert a feasible solution (feasible in terms of linear constraints only) into another feasible solution. At every iteration the algorithm considers active constraints only, the pressure on infeasible solutions is increased due to the decreasing values of temperature .

The method has an additional unique feature: it starts from a single point.2 Consequently, it is relatively easy to compare this method with other classical optimization methods whose performance are tested (for a given problem) from some starting point.

The method requires a starting and `freezing' temperatures, 0 and f , respectively, and the cooling scheme to decrease temperature . Standard values (reported by [START_REF] Michalewicz | Evolutionary optimization of constrained problems[END_REF]) are 0 = 1, i+1 = 0:1 i , with f = 0:000001.

One of the test problems [START_REF] Floudas | A Collection of Test Problems for Constrained Global Optimization Algorithms[END_REF] experimented by [START_REF] Michalewicz | Evolutionary optimization of constrained problems[END_REF] was:

minimize G6(x) = (x 1 10) 3 + (x 2 20) 3 , subject to nonlinear constraints: c1 :

(x 1 5) 2 + (x 2 5) 2 100 0, c2 :

(x 1 6) 2 (x 2 5) 2 + 82:81 0, and bounds: 13 x 1 100 and 0 x 2 100.

The known global solution is x = (14:095; 0:84296), and G6(x ) = 6961:81381 (see gure 3). Clearly, both constraints are active at the optimum. The starting point, which is not feasible, is x0 = (20:1; 5:84).

GENOCOP II approached the optimum very closely at the 12th iteration. The progress of the system is reported in the table 1.

In (Michalewicz, 1995a) some experimental evidence is presented to the eect that linear constraints of a problem may prevent the system to move closer to the optimum. This is an interesting example of damaging eect of limiting population to the feasible (with respect to linear constraints) region only. Additional experiments indicated that the method is very sensitive to the cooling scheme.

Methods of adaptive penalties

In this section we discuss two methods which are based on adaptive penalty functions.

The rst method was developed by [START_REF] Bean | A dual genetic algorithm for bounded integer programs[END_REF], and Hadj-Alouane and Bean (1992). As the previous method, it uses a penalty function, however, one component of the penalty function takes a feedback from the search process. Each individual is evaluated by the formula: feasible space optimum (14.095, 0.84296) Figure 3: A feasible space for test case of G6 eval(x) = f(x) + (t) P m j=1 f 2 j (x), where (t) is updated every generation t in the following way:

(t + 1) = 8 > > < > > :

(1= 1 ) (t); if bi 2 F for all t k + 1 i t 2 (t); if bi 2 S F for all t k + 1 i t (t); otherwise;

where bi denotes the best individual, in terms of function eval, in generation i, 1 ; 2 > 1 and 1 6 = 2 (to avoid cycling). In other words, the method (1) decreases the penalty component (t + 1) for the generation t + 1, if all best individuals in the last k generations were feasible, and (2) increases penalties, if all best individuals in the last k generations were infeasible. If there are some feasible and infeasible individuals as best individuals in the last k generations, (t + 1) remains without change.

The method introduces three additional parameters, 1 , 2 , and the time horizon, k. It seems that there is an interesting analogy between this method and the method used in early evolution strategies to optimize the convergence rate, where a \1/5 success rule" was proposed [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technisher Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Schoenauer | A survey of evolution strategies[END_REF]:

The ratio ' of successful mutations to all mutations should be 1/5. Increase the variance of the mutation operator, if ' is greater than 1/5; otherwise, decrease it. The 1/5 success rule emerged as a conclusion of the process of optimizing convergence rates of two functions (the so-called corridor model and sphere model). The rule was applied every k '(k) is the success ratio of the mutation operator during the last k generations, and c i > 1, c d < 1 regulate the increase and decrease rates for the variance of the mutation. Schwefel in his experiments [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] used the following values: c d = 0:82, c i = 1:22 = 1=0:82.

The intuitive reason behind adaptation of penalties in Bean's and Hadj-Alouane's method is similar to the reason behind the 1/5 success rule: the increased eciency of the search. In evolution strategies, if successful, the search would continue in \larger" steps; if not, the steps would be shorter. In Bean's and Hadj-Alouane's method, if constraints do not pose a problem, the search would continue with decreased penalties, if not, the penalties would be increased. The presence of both feasible and infeasible individuals in the set of best individuals in the last k generations means that the current value of penalty component (t) is set right.

The adaptive penalty function used by [START_REF] Smith | Genetic optimization using a penalty function[END_REF] incorporates both the search length and constraint severity feedback. It involves the estimation of a near-feasible threshold q j for each constraint 1 j m); such thresholds indicate distances from the feasible region F which are \reasonable" (or, in other words, which determine \interesting" infeasible solutions, i.e., solutions relatively close to the feasible region). Thus the evaluation function is dened as eval(x; t) = f(x) + F feas (t) F all (t) m X j=1 (f j (x)=q j (t)) k ;

where F all (t) denotes the unpenalized value of the best solution yet found (up to generation t), F feas (t) denotes the value of the best feasible solution yet found (up to generation t), and k is a constant. Note, that the near-feasible thresholds q j (t) are dynamic, i.e., they are adjusted during the search (for example, it is possible to dene q j (t) = q j (0)=(1 + j t) thus resulting in increasing the penalty component over time).

To the best of our knowledge, neither of the adaptive methods described in this subsection has been applied to continuous nonlinear programming problems.

Death penalty method

This method just rejects infeasible individuals (death penalty); the method has been used by evolution strategies (B ack et al., 1991) and simulated annealing. For some problems, this simple method provides quality results.

For example, one of the test cases considered recently (Michalewicz, 1995a) included the following problem [START_REF] Hock | Test Examples for Nonlinear Programming Codes[END_REF]. Minimize a function: G7(x) = x 2 1 + x 2 2 + x 1 x 2 14x 1 16x 2 + (x 3 10) 2 + 4(x 4 5) 2 + (x 5 3) 2 + 2(x 6 1) 2 + 5x 2 7 + 7(x 8 11) 2 + 2(x 9 10) 2 + (x 10 7) 2 + 45, subject to the following constraints: 105 4x 1 5x 2 + 3x 7 9x 8 0, 3(x 1 2) 2 4(x 2 3) 2 2x 2 3 + 7x 4 + 120 0, 10x 1 + 8x 2 + 17x 7 2x 8 0,

x 2 1 2(x 2 2) 2 + 2x 1 x 2 14x 5 + 6x 6 0, 8x 1 2x 2 5x 9 + 2x 10 + 12 0, 5x 2 1 8x 2 (x 3 6) 2 + 2x 4 + 40 0, 3x 1 6x 2 12(x 9 8) 2 + 7x 10 0, 0:5(x 1 8) 2 2(x 2 4) 2 3x 2 5 + x 6 + 30 0, and bounds 10:0 x i 10:0, i = 1; : : :; 10.

The problem has 3 linear and 5 nonlinear constraints; the function G7 where G7(x ) = 24:3062091. Six (out of eight) constraints are active at the global optimum (all except the last two).

The \death penalty" method provided with respectable solutions the best of which had the value of 25.653. To evaluate this method it is necessary to initialize a population by feasible solutions. This dierent initialization scheme makes the comparison of all the methods even harder. However, in experiments reported in (Michalewicz, 1995a), an interesting pattern emerged: the method generally gave a quite poor performance. Also, the method was not as stable as other methods; the standard deviation of returned solutions was relatively high.

Segregated genetic algorithm

The method of so-called segregated genetic algorithm was proposed by Le Riche et al. (1995) as yet another way to handle the problem of the robustness of the penalty level (see section 3.2). The tuning of the penalty level encounters the following dilemma: a too small penalty level leads to solutions which are infeasible (some penalized points still exhibit higher penalized tness than the best feasible point); a too high penalty level restricts the search inside the feasible region, forbidding any shot-cut across the infeasible region, and thus eventually failing to converge to the optimal solution.

The idea is to design two dierent penalized tness functions with static penalty terms p 1 and p 2 (see section 3.2): penalty p 1 is purposely too small, while penalty p 2 is hopefully too high. All individuals of the current population undergo crossover and mutation. The values of the two tness functions f i (x) = f(x) + p i (x), i = 1; 2, are computed for each resulting ospring (at no extra cost in term of objective function evaluation), and two ranked lists are created according to the values of the tnesses of all individuals (parents plus ospring) for each one of the tnesses. The selection of the parents of next generation is then achieved by picking up alternatively the best individual from each list (and removing this individual from both lists afterwards). The main idea is that such a selection scheme will result roughly in maintaining two subpopulations: the individuals selected on the basis of f 1 will more likely lie in the infeasible region while the ones selected on the basis of f 2 will probably stay in the feasible region; the overall process is thus allowed to reach the feasible optimum from both sides of the boundary of the feasible region.

This method gave excellent results in the domain of Laminated Design Optimization (Le Riche et al., 1995), but has not yet, to the best of our knowledge, been applied to continuous nonlinear programming problems.

Methods based on a search for feasible solutions

There are a few methods which emphasize the distinction between feasible and infeasible solutions in the search space S. One method considers the problem constraints in a sequence; a switch from one constraint to another is made upon arrival of a sucient number of feasible individuals in the population. The second method is based on an assumption, that any feasible solution is better than any infeasible one. The third method repairs infeasible individuals. We discuss these methods in turn.

Behavioral memory method

The method was proposed by [START_REF] Schoenauer | Constrained GA optimization[END_REF]; the authors called this method a \behavioral memory" approach. It works as follows: start with a random population of individuals (feasible or infeasible), set j = 1 (j is a constraint counter), evolve this population with eval(x) = f j (x), until a given percentage of the population (so-called ip threshold ) is feasible for this constraint,3 set j = j + 1, the current population is the starting point for the next phase of the evolution, where eval(x) = f j (x) (dened in the section 3.2). During this phase, points that do not satisfy one of the 1st, 2nd, ..., or (j 1)-th constraint are eliminated from the population. The stop criterion is again the satisfaction of the j-th constraint by the ip threshold percentage of the population. if j < m, repeat the last two steps, otherwise (j = m) optimize the objective function, i.e., eval(x) = f(x), rejecting infeasible individuals.

The method requires a linear order of all constraints which are processed in turn. It is unclear what is the inuence of the order of constraints on the results of the algorithm; our experiments indicated that dierent orders provide dierent results (dierent in the sense of the total running time and precision).

In total, the method requires three parameters: the sharing factor , the ip threshold , and a particular order of constraints. The method is very dierent to the other methods, and, in general, is dierent than other penalty approaches, since it considers only one constraint at the time. Also, in the last step of the algorithm the method optimizes the objective function f itself with a death penalty.

The method was evaluated by [START_REF] Schoenauer | Constrained GA optimization[END_REF] on the following problem: Maximize the function G8(x) = sin 3 (2x 1 )sin(2x 2 ) x 3 (x 1 +x 2 ) , subject to the following constraints: c1(x) = x 2 1 x 2 + 1 0, c2(x) = 1 x 1 + (x 2 4) 2 0 and bounds: 0 x 1 10 and 0 x 2 10. Figure 4: Function G8 together with two tness landscapes corresponding to two dierent penalization parameters : a small value of leaves the highest peaks outside the feasible region, while a large value of makes the feasible region look like a needle in a haystack (all plots are truncated between 10 and 10 for the sake of visibility) Function G8 has many local optima, the highest peaks are located along the x axis (e.g., G8(0:00015; 0:0225) > 1540). In the feasible region, however, G8 has two maxima of almost equal tness of value of 0.1. Figures 4 show the penalized tness landscape for dierent penalty parameters (i.e., plots of the function G8(x) (c1(x) + + c1(x) + ) for dierent values of ). All values of the penalized function are truncated at 10 and 10 to allow values in the feasible region (around 0) to be distinguishable.

This situation makes any attempt to use penalty parameters a dicult task: small penalty parameters leave the feasible region hidden among much higher peaks of the penalized tness (Figure 4-b) while penalty parameters large enough to allow the feasible region to really outstand in the penalized-tness landscape imply vertical slopes of that region (Figure 4-c): the partial credit principle, allowing the evolutionary algorithm to gently climb such slopes is not valid, and the discovery of the global maximum of the penalized tness thus relies on a lucky move from the huge and at low-tness region.

On the other hand, the behavioral memory method has no diculty in localizing 80% of the population rst into the feasible region for constraint c1 (step 1), then into the feasible region for both constraints (step 2), as can be seen on Figure 5. The optimization of G8 itself in that region is then straightforward (see the smooth bimodal tness landscape in that In other words, infeasible individuals have increased penalties: their values cannot be better than the value of the worst feasible individual (i.e., max x2F ff(x)g). 4For example, the above method was experimented with by Michalewicz (1995a) on one of the problems described by [START_REF] Hock | Test Examples for Nonlinear Programming Codes[END_REF]: minimize the function G9(x) = (x 1 10) 2 + 5(x 2 12) 2 + x 4 3 + 3(x 4 11) 2 + 10x 6 5 + 7x 2 6 + x 4 7 4x 6 x 7 10x 6 8x 7 , subject to the following constraint:

127 2x 2 1 3x 4 2 x 3 4x 2 4 5x 5 0, 282 7x 1 3x 2 10x 2 3 x 4 + x 5 0, 196 23x 1 x 2 2 6x 2 6 + 8x 7 0, 4x 2 1 x 2 2 + 3x 1 x 2 2x 2 3 5x 6 + 11x 7 0
and bounds: 10:0 x i 10:0, i = 1; : : : ; 7.

The problem has 4 nonlinear constraints; the function G9 is nonlinear and has its global minimum at

x = (2:330499; 1:951372; 0:4775414; 4:365726; 0:6244870; 1:038131; 1:594227),

where G9(x ) = 680:6300573. Two (out of four) constraints are active at the global optimum (the rst and the last one). The method of Powell and Skolnick gave a reasonable performance; the best solution found was 680.934. However, for some test problems (Michalewicz, 1995a) the method may have some diculties in locating a feasible solution.

Repairing infeasible individuals

This recently developed method, called Genocop III [START_REF] Michalewicz | Genocop III: A co-evolutionary algorithm for numerical optimization problems with nonlinear constraints[END_REF] is based on the idea of repairing infeasible solutions and also incorporates some concepts of co-evolution.

As in the original Genocop (section 3.1.1), linear equations are eliminated, number of variables is reduced, and linear inequalities are modied accordingly. All points included in the initial population satisfy linear constraints; specialized operators maintain their feasibility (in the sense of linear constraints) from one generation to the next. We denote the set of points which satisfy the linear constraints by F l S.

Nonlinear equations require an additional parameter () to dene the precision of the system. All nonlinear equations h j (x) = 0 (for j = q + 1; : : : ; m) are replaced by a pair of inequalities: h j (x) , so we deal only with nonlinear inequalities. These nonlinear inequalities restrict further the set F l : they dene the fully feasible part F F l of the search space S.

Genocop III incorporates the original Genocop system, but also extends it by maintaining two separate populations, where a development in one population inuences evaluations of individuals in the other population. The rst population P s consists of so-called search points from F l which satisfy linear constraints of the problem. As mentioned earlier, the feasibility (in the sense of linear constraints) of these points is maintained by specialized operators (see section 3.1.1). The second population P r consists of so-called reference points from F; these points are fully feasible, i.e., they satisfy all constraints.5 Reference points r from P r , being feasible, are evaluated directly by the objective function (i.e., eval(r) = f(r)). On the other hand, search points from P s are \repaired" for evaluation and the repair process works as follows. Assume, there is a search point s 2 P s . If s 2 F, then eval(s) = f(s), since s is fully feasible. Otherwise (i.e., s 6 2 F), the system selects6 one of the reference points, say r from P r and creates a sequence of random points z from a segment between s and r by generating random numbers a from the range h0; 1i: z = as + (1 a)r. 7 Once a fully feasible z is found, eval(s) = eval(z) = f(z). 8Additionally, if f(z) is better than f(r), then the point z replaces r as a new reference point in the population of reference points P r . Also, z replaces s in population of search points P s with some probability of replacement p r .

The structure of Genocop III is shown in gure 6 and the procedure for evaluating (not necessarily fully feasible) search points from population P s is given in gure 7. Note that there is some asymmetry between processing population of search points P s and population of reference points P r : while we apply selection procedure and operators to P s every generation, population P r is modied every k (parameter of the method) generations (however, some additional changes in P r are possible during the evaluation of search points, see gure 7). The main reason behind this arrangement is eciency of the system: search within feasible part of the search space F, as much less eective for NLP problems, is treated as a background event. Note also, that the \selection" and \alternation" steps are reversed in the evolution loop for the P r : due to a low probability of generating feasible ospring, rst parent individual reproduce and later the best feasible individuals (from both parents and ospring) are selected for the next population (this is similar to the ( + ) selection of evolution strategies (B ack et al., 1991).

procedure Genocop III begin t 0 initialize P s (t) initialize P r (t) evaluate P s (t) evaluate P r (t) while (not termination-condition) do begin t t + 1 select P s (t) from P s (t 1) alter P s (t) evaluate P s (t) if t mod k = 0 then begin alter P r (t) select P r (t) from P r (t 1) evaluate P r (t) end end end Figure 6: The structure of Genocop III Genocop III uses the objective function for evaluation of fully feasible individuals only, so the evaluation function is not distorted as in methods based on penalty functions. It introduces only few additional parameters (the population size of reference points, probability of replacement, frequency of application of operators to the population of reference points, precision ). It always returns a feasible solution. A feasible search space F is searched (population P r ) by making references from the search points and by application of operators (every some number of generations, k in gure 6). The neighborhoods of better reference points are explored more often. Some fully feasible points are moved into the population of search-points P s (replacement process), where they undergo additional transformation by specialized operators.

One of the most interesting parameters of the developed system is the probability of replacement p r (replacement of s by z in population of search points P s ; see gure 7). Recently, [START_REF] Orvosh | Shall we repair? Genetic algorithms, combinatorial optimization, and feasibility constraints[END_REF] reported a so-called 5%-rule: this heuristic rule states that in many combinatorial optimization problems, an evolutionary computation technique with a repair algorithm provides the best results when 5% of repaired individuals replace their infeasible originals. However, neither [START_REF] Davis | [END_REF] nor the authors are aware of any experiments connected with the probability of replacement for repair algorithms applied for problems in numerical domains. Experiments reported by [START_REF] Michalewicz | Genocop III: A co-evolutionary algorithm for numerical optimization problems with nonlinear constraints[END_REF] were done for p r = 0:20; this value procedure evaluate P s (t) begin for each s 2 P s (t) do if s 2 F then evaluate s (as f(s)) else begin select r 2 P r (t) generate z 2 F evaluate s (as f(z))

if f(r) > f(z) then replace r by z in P r replace s by z in P s with probability p r end end The results of experiments were quite good. For example, the following problem [START_REF] Hock | Test Examples for Nonlinear Programming Codes[END_REF] proved very dicult to all methods described earlier. The problem is to minimize G10(x) = x 1 + x 2 + x 3 , subject to the following constraints: 1 0:0025(x 4 + x 6 ) 0, 1 0:0025(x 5 + x 7 x 4 ) 0, 1 0:01(x 8 x 5 ) 0, x 1 x 6 833:33252x 4 100x 1 + 83333:333 0, x 2 x 7 1250x 5 x 2 x 4 + 1250x 4 0, x 3 x 8 1250000 x 3 x 5 + 2500x 5 0, and bounds 100 x 1 10000, 1000 x i 10000, i = 2; 3, 10 x i 1000, i = 4; : : :; 8.

The problem has 3 linear and 3 nonlinear constraints; the function G10 is linear and has its global minimum at x = (579:3167; 1359:943; 5110:071; 182:0174; 295:5985; 217:9799; 286:4162; 395:5979), where G10(x ) = 7049:330923. All six constraints are active at the global optimum.

For the above problem the best result of Genocop III was 7286.650: much better than the best result of the best system from those discussed in (Michalewicz, 1995a). Similar performance was observed on two other problems, G9 (with 680.640) and G7 (with 25.883). Additional interesting observation was connected with stability of the system. Genocop III had a very low standard deviation of results. For example, for problem G9, all results were between 680.640 and 680.889; on the other hand, other systems produced variety of results (between 680.642 and 689.660, see Michalewicz (1995a)).

Of course, all resulting points x were feasible, which was not the case with other systems (e.g., Genocop II produced a value of 18.917 for the problem G7, the systems based on the methods of Homaifar, Lai, and Qi (section 3.2.1) and Powell and Skolnick (section 3.3.2) gave results of 2282.723 and 2101.367, respectively, for the problem G10, and these solutions were not feasible).

Hybrid methods

It is relatively easy to develop hybrid methods which combine evolutionary computation techniques with deterministic procedures for numerical optimization problems. [START_REF] Waagen | The stochastic direction set algorithm: A hybrid technique for nding function extrema[END_REF] combined evolutionary programming technique (with oating point representation, Fogel, 1995) with the direction set method of Hooke-Jeeves; this hybrid method was tested on three (unconstrained) test functions. [START_REF] Myung | Preliminary investigation into a two-stage method of evolutionary optimization on constrained problems[END_REF] considered a similar approach, but they experimented with constrained problems. Again, they combined oating-point evolutionary programming technique with some other method|developed by [START_REF] Maa | A two-phase optimization neural network[END_REF]. However, while the method of [START_REF] Waagen | The stochastic direction set algorithm: A hybrid technique for nding function extrema[END_REF] incorporated the direction set algorithm as a problem-specic operator of his evolutionary technique, [START_REF] Myung | Preliminary investigation into a two-stage method of evolutionary optimization on constrained problems[END_REF] divided the whole optimization process into two separate phases. During the rst phase, evolutionary algorithm optimizes the function eval

(x) = f(x) + s 2 0 @ m X j=1 f 2 j (x) 1 A ;
where s is a constant. After the termination of this phase, the second phase of the optimization algorithm of [START_REF] Maa | A two-phase optimization neural network[END_REF] was applied to the best solution found during the rst phase; this phase iterates until the system

x0 = 5 f(x) 2 4 m X j=1 5f j (x)(sf j (x) + j ) 3 5 ;
is in equilibrium, where the Lagrange multipliers are updated as 0 j = sf j for a small positive constant . This hybrid method was successfully applied to a few test cases; for example, one of the test cases was to minimize G11(x) = x 2 1 + (x 2 1) 2 , subject to a nonlinear constraint:

x 2 x 2 1 = 0, and bounds:

1 x i 1, i = 1; 2.
The known global solutions are x = (0:70711; 0:5), and G11(x ) = 0:75000455. For this problem, the rst phase of the evolutionary process converges quickly in 100 generations to (0:70711; 0:5), and the second phase drives the system to the global solution. However, the method was not tested for problems of higher dimensions.

Several other constraint handling methods deserve also some attention. For example, some methods use of the values of objective function f and penalties f j (j = 1; : : : ; m) as elements of a vector and apply multi-objective techniques to minimize all components of the vector. For example, [START_REF] Schaer | Multiple objective optimization with vector evaluated genetic algorithms[END_REF] Vector Evaluated Genetic Algorithm (VEGA) selects 1=(m + 1) of the population based on each of the objectives. Such an approach was incorporated by [START_REF] Parmee | The development of directed genetic search technique for heavily constrained design spaces[END_REF] in the development of techniques for constrained design spaces. On the other hand, in the approach by [START_REF] Surry | A multi-objective approach to constrained optimization of gas supply networks[END_REF], all members of the population are ranked on the basis of constraint violation. Such rank r, together with the value of the objective function f, leads to the two-objective optimization problem. This approach gave a good performance on optimization of gas supply networks [START_REF] Surry | A multi-objective approach to constrained optimization of gas supply networks[END_REF].

Also, an interesting approach was reported by [START_REF] Paredis | Coevolutionary constraint satisfaction[END_REF]. The method (described in the context of constraint satisfaction problems) is based on a co-evolutionary model, where a population of potential solutions co-evolves with a population of constraints: tter solutions satisfy more constraints, whereas tter constraints are violated by more solutions. There is some development connected with generalizing the concept of \ant colonies" [START_REF] Colorni | Distributed optimization by ant colonies[END_REF] (which were originally proposed for order-based problems) to numerical domains [START_REF] Bilchev | Ant colony search vs. genetic algorithms[END_REF]; rst experiments on some test problems (from Table 3) gave very good results [START_REF] Bilchev | [END_REF]. It is also possible to incorporate the knowledge of the constraints of the problem into the belief space of cultural algorithms [START_REF] Reynolds | An introduction to cultural algorithms[END_REF]; such algorithms provide a possibility of conducting an ecient search of the feasible search space [START_REF] Reynolds | Using cultural algorithms for constraint handling in Genocop[END_REF]. Recently, [START_REF] Kelly | [END_REF] developed a new procedure for dicult constrained optimization problems, which \...has succeeded in generating optimal and near-optimal solutions in minutes to complex nonlinear problems that otherwise have required days of computer time in order to verify optimal solutions. These problems include nonconvex, nondierentiable and mixed discrete functions over constrained regions. The feasible space explicitly incorporates linear and integer programming constraints and implicitly incorporates other nonlinear constraints by penalty functions. The objective functions to be optimized may have no `mathematically speciable' representation, and can even require simulations to generate. Due to this exibility, the approach applies also to optimization under uncertainty (which increases its relevance to practical applications)."

The method combines evolutionary processes with (1) scatter search, (2) with procedures of mixed integer optimization, and (3) with adaptive memory strategies of tabu search.

However, all the above methods are on early stages of their development, or they have not been applied yet to the general NLP, so we do not discuss them any further.

Conclusions and Future Work

It is dicult to compare the constraint-handling methods which were discussed in the previous section of this paper. The main reason is that the original formulations of these methods are based on dierent representations (binary or oat), and consequently they use dierent operators. For example, the method of static penalties (section 3.2.1) was proposed for a standard GA with binary representation, standard crossover and mutation. Genocop (section 3.1.1) uses oating point representation with seven specialized operators. The method of superiority of feasible points (section 3.3.2) was also developed for oating point representation. However, their operators were dierent to those of Genocop. Bean and Hadj-Alouane adapted penalties (section 3.2.4) for integer programming problems and used binary representation to code integers. Death penalty methods (section 3.2.5) were mainly used by evolution strategies (oating point representation with Gaussian mutation and adaptation of control parameters), and the segregated GAs (section 3.2.6) has not been applied yet to the NLP. For these reasons we do not provide numerical results of the above methods on all test cases dened in the previous section. The only general observation we can make at this stage is that for constrained numerical optimization problems, the oating point representation provides with much better results (some of experiments on comparisons between binary and oating point representations for unconstrained problems were reported in Michalewicz, 1996, chapter 5).

Some partial comparisons of several constraint-handling methods were reported (Michalewicz, 1995a), where a modied version of Genocop system was used for all experiments; the code was updated accordingly for each considered method. The results of these comparisons are given in Table 2.

It is not clear what characteristics of a constrained problem make it dicult for an evolutionary technique. It seems that no single parameter (number of linear, nonlinear, active constraints, the ratio = jFj=jSj of feasible points in S, type of the function, number of variables) proved to be signicant as a major measure of diculty of the problem. For example, many methods approached the optimum quite closely for the test cases G1 and G7 (with = 0:0111% and = 0:0003%, respectively), whereas most of the methods experienced diculties for the test case G10 (with = 0:0010%). Two quadratic functions (the test cases G1 c 0, 1, 1 0, 0, 0 0, 1, 0 1, 0, 0 0, 0, 0 Table 2: Experimental results. For some methods, the best (b), median (m), and the worst (w) result (out of 10 independent runs) are reported together with the number (c) of violated constraints at the median solution: the sequence of three numbers indicate the number of violations by more than 1.0, more than 0.1, and more than 0.001, respectively. The symbol `|' stands for `the results were not meaningful'. The results for method of section 3.2.5 report also the results of experiments (case 3.2.5(f)) where the initial population was feasible.

and G7) with a similar number of constraints (9 and 8, respectively) and an identical number (6) of active constraints at the optimum, gave a dierent challenge to most of these methods. Also, several methods were quite sensitive to the presence of a feasible solution in the initial population (e.g., method of section 3.2.5; see table 2). All test cases are summarized in table 3; for each test case we list number n of variables, type of the function f, the relative size of the feasible region in the search space given by the ratio , the number of constraints of each category (linear inequalities LI, nonlinear equations NE and inequalities NI), and the number a of active constraints at the optimum (including equality constraints).

Despite some obvious weaknesses of the method of static penalties (section 3.2.1), this basic method, which requires xed user-supplied penalty parameters, remains the most popular one to handle constraints for evolutionary algorithms and many other optimization frameworks. The main reason for that popularity is that it certainly is the simplest technique to implement: it requires only a straightforward modication of the objective function. Moreover, it is fair to say that it gives quite reasonable results in many cases. Nevertheless, as emphasized in some of Function n Type of f LI NE NI a G1 13 quadratic 0.0111% 9 0 0 6 G2 k nonlinear 99.8474% 0 0 2 1 G3 k polynomial 0.0000% 0 1 0 1 G4 5 quadratic 52.1230% 0 0 6 2 G5 4 cubic 0.0000% 2 3 0 3 G6 2 cubic 0.0066% 0 0 2 2 G7 10 quadratic 0.0003% 3 0 5 6 G8 2 nonlinear 0.8560% 0 0 2 0 G9 7 polynomial 0.5121% 0 0 4 2 G10 8 linear 0.0010% 3 0 3 6 G11 2 quadratic 0.0000% 0 1 0 1 Table 3: Summary of eleven test cases. The ratio = jFj=jSj was determined experimentally by generating 1,000,000 random points from S and checking whether they belong to F (for G2

and G3 we assumed k = 50). LI, NE, and NI represent the number of linear inequalities, and nonlinear equations and inequalities, respectively the research cited in this paper, there are also many other cases for which the penalty method fails. This raises two questions:

(1) what features of the problem can guide the choice of the constraint-handling technique to be incorporated into an evolutionary algorithm (either among those surveyed here, or based on a new paradigm), and (2) is it possible to detect a priori such problems, which are resistant to the penalty approach? Let us provide a few comments addressing these two questions.

The experimental results of this paper suggest that making an appropriate a priori choice of an evolutionary method for a nonlinear optimization problem remains an open question. It seems that more complex properties of the problem (e.g., the characteristic of the objective function together with the topology of the feasible region) may constitute quite signicant measures of the diculty of the problem. Also, some additional measures of the problem characteristics due to the constraints might be helpful. However, this kind of information is not generally available a priori. One possibility would be to measure how fast the feasible region grows when the constraints are relaxed. We discuss briey these concepts in the following paragraphs.

Consider the general NLP problem of section 1. Suppose that each constraint g j ; j 2 [1; q] and h j ; j 2 [q + 1; m] is upper-bounded on S by M j 0. Then, for every value of ! 2 [0; 1], we can dene F ! as F ! = fx 2 S : 8j 2 [1; q] g j (x) !M j and 8j 2 [q + 1; m] jh j (x)j !M j g: Of course, the diculty of the NLP also heavily depends on the objective function f itself. So one possible measure of diculty of the constrained problem is an evaluation of the optimum of the objective function f with relaxed constraints. The following formula,9 ' ! = max x2F f(x) max x2F f(x) ; can serve as a such measure. Unfortunately, this ratio can only be estimated a posteriori by solving a large number of NLPs, all as dicult as the original one.

In general, there are three possible sources of diculties of the NLP. These diculties can arise from the objective function itself, from the denition of the feasible region (i.e., from the constraints alone), or from the coupling between the objective function and the constraints. Hence the following features should be taken into account when trying to characterize dierent NLPs:

The ruggedness of the unconstrained tness landscape certainly is the rst important characteristic inuencing the overall problem diculty; for instance, linear or convex objective functions will always result in easier constrained problems than almost-chaotic functions, for the same set of constraints; however, this characteristic can hardly be estimated a priori.

The sparseness of the feasible region indeed is a crucial factor of diculty: in some problems, nding a single feasible point is the major diculty. Moreover, a high slope of the constraints on the border of the feasible region is another way by which a constrained problem can resist the penalty method; some hints have been given above (the factors and ! can be numerically estimated a priori).

A high ratio between the highest global optima of the objective function (on the whole domain where it is dened) and of the optima of the constrained function (e.g., the local optima of the objective function in the feasible region), as well as a small distance between these global optima and the feasible region can make the constrained optimization problem almost intractable for penalty methods; unfortunately, it is almost impossible to estimate this factor a priori. However, some useful a posteriori informations about the behavior of dierent methods might be obtained at the cost of heavy computations by comparing their respective performances on gradually relaxed hard NLPs.

The number of active constraints at the optimum is of course of importance: the more constraints that are active at the optimum, the more likely to succeed are algorithms searching close to the boundary of the feasible region; the extreme case being that described in section 3.1.2. But here again it is very dicult to predict such a number a priori.

It seems that the most promising approach at this stage of research is experimental, involving the design of a scalable test suite of constrained optimization problems, in which many of these features could be easily tuned. Then it should be possible to test new methods with respect to the corpus of all available methods.

Such experimental suite should contain:

Convex problems (for validation only, as any gradient-bassed method will outperform any stochastic algorithm by several orders of magnitude on such problems);

Convex objective function and non-convex constraints (from easy to dicult in slope, number of connected components, etc.);

Non-convex objective function (with multiple optima) with convex constraints; Non-convex objective function (with numerous multiple optima) and highly non-convex steep constraints (with F consisting of possibly disjoint regions).

Such an experimental test suite is currently under construction.

Figure 1 :

 1 Figure 1: The graph of function G2 for n = 2. Infeasible solutions were assigned value zero

  Figure2: The geometrical crossover on the hyperboloid around its center: when parameter goes from 0 to 1, the ospring of two parents denes the line joining them on the plot.

Figure 7 :

 7 Figure 7: Evaluation of population P s

  It is interesting to observe that F 0 = F and F 1 = S. The way ! = jF ! j=jSj behaves around ! = 0 | measured, for instance, by the value! = d ! d! j !=0 ;can serve as an indicator of \how fast" the feasible region expands around its boundary.

Table 1 :

 1 Progress of GENOCOP II on test case of function G6; for iteration 0 the best point is the starting point.

		Iteration number 0 1 2 3 4 5 6 7 8 9 10 11 12	The best point (20.1, 5.84) (13.0, 0.0) (13.63, 0.0) (13.63, 0.0) (13.73, 0.16) (13.92, 0.50) (14.05, 0.75) (14.05, 0.76) (14.05, 0.76) (14.10, 0.87) (14.10, 0.86) (14.10, 0.85) (14.098, 0.849)	Active constraints c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2 c 1 , c 2
	generations:	t+1 =	8 > < > :	c d t ; if '(k) < 1=5; c i t ; if '(k) > 1=5; t ; if p s (k) = 1=5;
	where

This feature, however, is not essential. The only important requirement is that the next population contains the best individual from the previous population.

The method suggests the use of a sharing scheme to maintain diversity of the population.

[START_REF] Powell | Using genetic algorithms in engineering design optimization with non-linear constraints[END_REF] achieved the same result by mapping evaluations of feasible solutions into the interval ( 1; 1) and infeasible solutions into the interval (1; 1). For ranking and tournament selections this implementational dierence is not important.

If Genocop III has diculties in locating such a reference point for the purpose of initialization, the user is prompted for it. In cases, where the ratio jFj=jSj of feasible points in the search space is very small, it may happen that the initial set of reference points consists of a multiple copies of a single feasible point.

Better reference points have better chances to be selected; a selection method based on nonlinear ranking was used.

Note that all such generated points z belong to F l .

Clearly, in dierent generations the same search point S can evaluate to dierent values due to the random nature of the repair process.

For maximization problems. The objective function f takes non-negative values.

Acknowledgments:

This material is based upon work supported by the National Science Foundation under Grant IRI-9322400. The authors would to thank Ken De Jong for his suggestions which made the article more readable.

Figure 5: Illustration of a sample run of behavioral memory method on G8: zoom around the feasible region. The dashed curve is constraint c1 and the continuous curve is constraint c2.

The population at end of rst step samples the feasible region for c1 while that of the end of step 2 samples the feasible region for both c1 and c2 region in Figure 4). Throughout all steps, however, the diversity in the population had to be preserved: the resulting population is to be used as the starting point for another evolution, and hence must sample as uniformly as possible the target feasible region in the two rst steps; the objective function is multi-modal in that feasible region, hence care must be taken to locate both peaks. This requirement, achieved using the sharing scheme, is even more mandatory when the feasible region is made of several distinct connected components, as shown by [START_REF] Schoenauer | Constrained GA optimization[END_REF] on a similar problem (with a slightly dierent constraint c2).

Method of superiority of feasible points

The method was developed by [START_REF] Powell | Using genetic algorithms in engineering design optimization with non-linear constraints[END_REF]. The method is based on a classical penalty approach with one notable exception. Each individual is evaluated by the formula: eval(x) = f(x) + r P m j=1 f j (x) + (t;x), where r is a constant; however, the original component (t;x) is an additional iteration dependent function which inuences the evaluations of infeasible solutions. The point is that the method distinguishes between feasible and infeasible individuals by adopting an additional heuristic rule (suggested earlier by [START_REF] Richardson | Some guidelines for genetic algorithms with penalty functions[END_REF]: for any feasible individual x and any infeasible individual ỹ: eval(x) < eval(ỹ), i.e., any feasible solution is better than any infeasible one. This can be achieved in many ways; one possibility is to set (t;x) = ( 0; if x 2 F maxf0; max x2F ff(x)g min x2S F ff(x) + r P m j=1 f j (x)gg; otherwise:


Ivana Kojčić
email: ivana.kojcic@inria.fr

Théodore Papadopoulo

Rachid Deriche

Samuel Deslauriers-Gauthier

Connectivity

Connectivity-informed M/EEG inverse problem

Keywords: MEG, source localization, inverse problems, diffusion MRI

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Information between brain regions is transferred through white matter fibers with delays that are measurable with magnetoencephalography and electroencephalography (M/EEG) due to its millisecond temporal resolution. Therefore, a useful representation of the brain is that of a graph where its nodes are the cortical areas and edges are the physical connections between them: either local (between adjacent vertices on the cortical mesh) or non-local (long-range white matter fibers). These long-range anatomical connections can be obtained by diffusion MRI (dMRI) tractography, thus giving us an insight on interaction delays of the macroscopic brain network. A fundamental role in shaping the rich temporal structure of functional connectivity is played by the structural connectivity [START_REF] Sporns | The human connectome: a complex network[END_REF] that places constraints on which functional interactions occur in the network. In the context of regularizing the dynamics of M/EEG and recovering electrical activity of the brain from M/EEG measurements, traditional linear inverse methods deploy different constraints such as minimum norm, maximum-smoothness in space and/or time along the cortical surface. However, they usually do not take into account the structural connectivity and very few include delays supported by dMRI as a prior information [START_REF] Deslauriers-Gauthier | White matter information flow mapping from diffusion mri and eeg[END_REF]. The goal of this work is to include these delays into the MEG source reconstruction process by imposing temporal smoothness in structurally connected sources, with the corresponding delays. We propose to encapsulate delays provided by dMRI in a graph representation and show their potential in improving the MEG source reconstruction when compared to a state-of-the-art approach [START_REF] Hammond | Cortical graph smoothing: a novel method for exploiting dwi-derived anatomical brain connectivity to improve eeg source estimation[END_REF].

Forward and inverse M/EEG problem

Distributed source models place the current sources at a large number of vertices distributed on the cortex. The relationship between source amplitudes and M/EEG measurements is expressed by the linear model M = GJ + E. The matrix of measurements is given by M ∈ IR N ×T with N sensors and T time samples. The unknown matrix of S source amplitudes is given by J ∈ IR S×T . The gain (lead field) matrix G ∈ IR N ×S provides a linear relationship between source amplitudes and sensor data (a.k.a the M/EEG forward solution) and E ∈ IR N ×T is additive noise. Recovering electrical activity of the brain from M/EEG data is ill-posed and prior assumptions need to be introduced. A wide range of generalizations of minimum norm estimators [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: minimum norm estimates[END_REF] exist in the form

min J U (J ) = min J {F (M, J ) + λP (J )} = min J M -GJ 2 2 + λ W J 2 2 (1)
which have a closed-form solution

Ĵ = G T G + λW T W -1 G T M . They differ
by the choice of the matrix W ∈ IR S×S which incorporates a certain regularization property. However, all l 2 -based inverse solvers suffer from smearing of even focal activation and often fail to exploit specific knowledge about the brain, such as the temporal dynamics of the underlying sources.

Laplacian as a smoothness constraint We can consider brain as a directed graph where its nodes are the current sources and its edges encode physical connection between them: either local (between adjacent vertices on the cortical mesh) or non-local (long-range white matter fibers). In the LORETA approach [START_REF] Pascual-Marqui | Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain[END_REF], similar activity is favoured between spatially adjacent vertices on the cortical mesh by choosing W to be the Laplacian matrix on the cortical surface, giving the solution with the maximum spatial smoothness. Since this approach does not account for long-range connectivity, it was extended in the Cortical Graph Smoothing (CGS) method [START_REF] Hammond | Cortical graph smoothing: a novel method for exploiting dwi-derived anatomical brain connectivity to improve eeg source estimation[END_REF] by forming a hybrid local/nonlocal connectivity graph with A loc (spatial adjacency on the cortical surface) and A tr (tractography-based connectivity matrix). They are used to form local L loc and tractography-based L tr graph Laplacians, resulting in the penalty term P CGS (J ) = λ loc J T L loc J + λ tr J T L tr J . This way CGS penalizes the weighted sum of squared differences in activity between connected cortical patches.

Connectivty-Informed MEG Inverse Problem (CIMIP)

While previous approaches exploited structural connectivity between different cortical areas, transmission delays were not taken into account. In order to enforce temporal smoothness between the time courses of connected sources, we incorporated a modified version of the Laplacian operator as a penalty term in the minimization. We can consider all time samples in a single large problem as m = Gj + where m ∈ IR N T and j ∈ IR ST are vectors of concatenated measurements and source intensities, G = diag(G, ...G) ∈ IR N T ×ST and ∈ IR N T is additive noise. To recover electrical activity of the brain, we minimize the following objective function

U CIMIP (j) = m -Gj 2 2 + λ(j T Lj) 2 (2)
Firstly, a binary time-dependent connectivity graph à = Ãloc + Ãtr is built. Short-range connections are encapsulated in Ãloc = diag(A loc , ...A loc ) ∈ IR ST ×ST . In order to include connections for different delays, we build a tractographybased graph Ãtr ∈ IR ST ×ST where nonzero elements designate the presence of a long-range connection for a specific delay, which was not considered in the CGS. The output of dMRI tractography is a set of streamlines representing the white matter fiber bundles. Given the streamline lengths and the information conduction speed, we can calculate the propagation delays for each connection.

In accordance to recent findings [START_REF] Drakesmith | Estimating axon conduction velocity in vivo from microstructural mri[END_REF] showing that the proportionality between fiber diameter and conduction velocity is 6.67ms -1/µm and previous works of [START_REF] Deslauriers-Gauthier | White matter information flow mapping from diffusion mri and eeg[END_REF], information conduction speed is assumed to be constant across brain and equal to 6m/s. Using Ã, the matrix L in ( 2) is a regularization matrix defined as L = I -Ã. This is a symmetric block matrix where each block corresponds to a specific delay. In the standard Laplacian the sum of each row is 0, so the preferred solutions would be the ones where all connections of a given node have the same amplitude. For that reason, we used L in order to favour solutions where the activity of a source is equal to the sum of activities of its connections.

Data simulation

Source-level brain activity is modelled as a function of its local and long-range connections as a multivariate autoregressive (MAR) model. The source space was parcellated into 68 neuroanatomical regions of interest. Assuming that only a few regions are typically active during a cognitive task, the activity is simulated for 10 randomly chosen subnetworks as a graph path with different lengths, so the number of active regions varied from 2 to 5. Two types of simulations were performed: multiple focal sources and multiple spread sources. In the first, only the sources in the active region that have a long-range connection in the next active region were assigned the activity. In the latter, this activity was also spread to its adjacent neighbours. The activity was propagated between sources that share a long-range connection, with an appropriate delay.

Results

Results of the reconstruction of brain activity from simulated MEG measurements were compared to the CGS method. Source estimation performance was evaluated with the same metric as in [START_REF] Hammond | Cortical graph smoothing: a novel method for exploiting dwi-derived anatomical brain connectivity to improve eeg source estimation[END_REF]

Conclusion

Our preliminary results show that including the conduction delays provided by dMRI encapsulated in a graph representation improves MEG source localization when compared to the CGS method. Future work will examine the performance on real MEG data and comparison to other methods.

Table 1 .

 1 , by identifying the largest local maxima of the estimated current sources. Peak localization error was computed as the Euclidean distance between local spatial maxima k * of simulated (GT) and reconstructed (REC) sources i.e. d(k * ) = |k * REC -k * GT | Mean and standard deviation of peak localization errors (in mm)

	Method Multiple focal sources Multiple spread sources
	CGS	48.09 ± 17.98	53.03 ± 17.2
	CIMIP 28.05 ± 8.4	27.99 ± 6.97

Acknowledgments -this work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 649665: CoBCoM -Computational Brain Connectivity Mapping)