A Multiagent Evolutionary Algorithm for Constraint ... - CiteSeerX

2 downloads 0 Views 3MB Size Report
are integrated to solve global numerical optimization problems, and the method can find high quality solutions at a low. A Multiagent Evolutionary Algorithm for.
SMCB-E-03172004-0141.R2

1

A Multiagent Evolutionary Algorithm for Constraint Satisfaction Problems Weicai Zhong, Jing Liu, and Licheng Jiao, Senior Member, IEEE

Abstract—With the intrinsic properties of constraint satisfaction problems (CSPs) in mind, we divide CSPs into two types, namely, permutation CSPs and non-permutation CSPs. According to their characteristics, several behaviors are designed for agents by making use of the ability of agents to sense and act on the environment. These behaviors are controlled by means of evolution, so that the multiagent evolutionary algorithm for constraint satisfaction problems (MAEA-CSPs) results. To overcome the disadvantages of the general encoding methods, the minimum conflict encoding is also proposed. Theoretical analyses show that MAEA-CSPs has a linear space complexity and converges to the global optimum. The first part of the experiments uses 250 benchmark binary CSPs and 79 graph coloring problems from the DIMACS challenge to test the performance of MAEA-CSPs for non-permutation CSPs. MAEA-CSPs is compared with six well-defined algorithms and the effect of the parameters is analyzed systematically. The second part of the experiments uses a classical CSP, n-queen problems, and a more practical case, job-shop scheduling problems (JSPs), to test the performance of MAEA-CSPs for permutation CSPs. The scalability of MAEA-CSPs along n for n-queen problems is studied with great care. The results show that MAEA-CSPs achieves good performance when n increases from 104 to 107, and has a linear time complexity. Even for 107-queen problems, MAEA-CSPs finds the solutions by only 150 seconds. For JSPs, 59 benchmark problems are used, and good performance is also obtained. Index Terms—Constraint satisfaction problems, evolutionary algorithms, graph coloring problems, job-shop scheduling problems, multiagent systems, n-queen problems.

A

I. INTRODUCTION

large number of problems coming from artificial intelligence as well as other areas of computer science and engineering can be stated as Constraint Satisfaction Problems (CSPs) [1]-[3]. A CSP has three components: variables, values, and constraints. The purpose is to find an assignment of values to variables such that all constraints are satisfied. Many methods have been proposed to deal with CSPs, where backtracking and generate-and-test methods are well-known traditional approaches [1]. Although the backtracking method eliminates a

Manuscript received March 17, 2004. This work is supported by the National Natural Science Foundation of China under Grant 60133010 and 60372045, and the “863” project under Grant 2002AA135080. The authors are with the Institute of Intelligent Information Processing, Xidian University, Xi’an, 710071, China Jing Liu, corresponding author, phone: 86-029-88202661; fax: 86-029-88201023; e-mail: [email protected].

subspace from the Cartesian product of all the variable domains, the computational complexity for solving most nontrivial problems remains exponential. Many studies have been conducted to investigate various ways of improving the backtracking method, such as consistency techniques [4], [5], the dependency-directed backtracking scheme [6], [7]. Nevertheless, it is still unable to solve nontrivial large-scale CSPs in a reasonable runtime. For the generate-and-test method, one of the most popular ideas is to employ a local search technique [8]-[10]. It generates an initial configuration and then incrementally uses “repair” or “hill climbing” to modify the inconsistent configuration to move to a neighborhood configuration having the best or better evaluation value among the neighbors, until a solution is found. As related to the idea of local search, other heuristics have also been developed, such as the min-conflicts heuristic [11], [12], GSAT [13]. In addition, local search techniques are often combined with Evolutionary Algorithms (EAs), another popular kind of method for solving CSPs. Historically, CSPs have been approached from many angles by EAs. In this field, the method used to treat the constraints is very important. The constraints can be handled directly, indirectly or by mixed methods [14]. Among the available methods, some put the emphasis on the usage of heuristics, such as ARC-GA [15], [16], COE-H GA [17], [18], Glass-Box [19], H-GA [20], [21], whereas others handle the constraints by fitness function adaptation, such as CCS [22], [23], MID [24]-[26], SAW [27], [28]. These methods dealt with CSPs by EAs from different angles and boosted the development of this field. Reference [29] has made an extensive performance comparison of all these EAs on a systematically generated test suite. Agent-based computation has been studied for several years in the field of distributed artificial intelligence [30], [31] and has been widely used in other branches of computer science [32]-[35]. Problem solving is an area with which many multiagent-based applications are concerned. It includes distributed solutions to problems, solving distributed problems, and distributed techniques for problem solving [30], [31]. Reference [33] introduced an application of distributed techniques for solving constraint satisfaction problems. They solved 7000-queen problems by an energy-based multiagent model. In [34], multiagent systems and genetic algorithms (GAs) are integrated to solve global numerical optimization problems, and the method can find high quality solutions at a low

SMCB-E-03172004-0141.R2 computational cost even for the functions with 10 000 dimensions. All these results show that both agents and EAs have a high potential in solving complex and ill-defined problems. In this paper, with the intrinsic properties of CSPs in mind, we divide CSPs into two types, namely, permutation CSPs and non-permutation CSPs. According to the different characteristics of the two types of problems, several behaviors are designed for agents. Furthermore, all such behaviors are controlled by means of evolution, so that a new algorithm, multiagent evolutionary algorithm for constraint satisfaction problems (MAEA-CSPs), results. In MAEA-CSPs, all agents live in a latticelike environment. Making use of the designed behaviors, MAEA-CSPs realizes the ability of agents to sense and act on the environment in which they live. During the process of interacting with the environment and the other agents, each agent increases energy as much as possible, so that MAEA-CSPs can find solutions. Experimental results show that MAEA-CSPs provides good performance. From another viewpoint, since MAEA-CSPs uses a lattice-based population, it can be also considered as a kind of fined-grained parallel GA. Fined-grained parallel GAs are a kind of parallel implementation of GAs [36]-[45], and have been applied to many problems, such as combinatorial optimization [46], function optimization [47], [48], scheduling problems [49]. But since they are just parallel techniques, they often present the same problem of premature convergence as traditional GAs [50]. MAEA-CSPs makes use of the ability of agents in sensing and acting on the environment, and puts emphasis on designing behaviors for agents. The experimental results show that MAEA-CSPs achieves good performance for various CSPs, that is, binary CSPs, graph coloring problems (GCPs), n-queen problems, and job-shop scheduling problems (JSPs). Therefore, MAEA-CSPs is a successful development of fined-grained parallel GAs. The rest of this paper is organized as follows: Section II describes the agents designed for CSPs. Section III describes the implementation of MAEA-CSPs, and analyzes the space complexity and convergence. Section IV uses binary CSPs and GCPs to investigate the performance of MAEA-CSPs on non-permutation CSPs, whereas Section V uses n-queen problems and JSPs to investigate the performance on permutation CSPs. Finally, conclusions are presented in Section VI. II. CONSTRAINT SATISFACTION AGENTS According to [31] and [33], an agent is a physical or virtual entity essentially having the following properties: (a) it is able to live and act in the environment; (b) it is able to sense the local environment; (c) it is driven by certain purposes and (d) it has some reactive behaviors. Multiagent systems are computational systems in which several agents interact or work together in order to achieve purposes. As can be seen, the meaning of an agent is very comprehensive, and what an agent represents is different for different problems. In general, four elements

2 should be defined when multiagent systems are used to solve problems. The first is the meaning and the purpose of each agent. The second is the environment in which all agents live. Since each agent has only local perceptivity, so the third is the local environment. The last is the behaviors that each agent can take to achieve the purpose. A. Constraint satisfaction problems A CSP has three components: 1) A finite set of variables, x={x1, x2, …, xn}; 2) A domain set D, containing a finite and discrete domain for each variable:

{

}

D={D1, D2, …, Dn}, xi∈Di= d1 , d 2 ,  , d| Di | , i=1, 2, …, n (1) where |•| stands for the number of elements in the set; 3) A constraint set, C={C1(x1), C2(x2), …, Cm(xm)}, where xi, i=1, 2, …, m is a subset of x, and Ci(xi) denotes the values that the variables in xi cannot take simultaneously. For example, given a constraint C({x1, x2})=〈d1, d2〉, it means when x1=d1, d2 cannot be assigned to x2, and when x2=d2, d1 cannot be assigned to x1. Thus, the search space of a CSP, S, is a Cartesian product of the n sets of finite domains, namely, S=D1×D2×…×Dn. A solution for a CSP, s=〈s1, s2, …, sn〉∈S, is an assignment to all variables such that the values satisfy all constraints. Here is a simple example: Example 1: A CSP is described as follows:  x = {x1 , x2 , x3 } ⌠ ⌠ D = {D1 , D2 , D3 }, Di = {1, 2, 3}, i = 1, 2,3 ⌠C = {C ({x , x }) = ∑1,3⌡, C ({x , x }) = ∑3,3⌡, 1 1 2 2 1 2 ⌠ ⌠ C3 ({x1 , x3 }) = ∑2,1⌡, C4 ({x1 , x3 }) = ∑2,3⌡, ⌠ (2) 〉 C5 ({x1 , x3 }) = ∑3,1⌡, C6 ({x1 , x3 }) = ∑3,3⌡, ⌠ ⌠ C7 ({x2 , x3 }) = ∑1,1⌡, C8 ({x2 , x3 }) = ∑1, 2⌡, ⌠ ⌠ C9 ({x2 , x3 }) = ∑1,3⌡, C10 ({x2 , x3 }) = ∑2,1⌡, ⌠ C11 ({x2 , x3 }) = ∑3,1⌡} ⌠∫ All solutions for this CSP are 〈1, 2, 2〉, 〈1, 2, 3〉, 〈2, 2, 2〉, 〈2, 3, 2〉, 〈3, 2, 2〉.  One may be interested in finding one solution, all solutions or in proving that no solution exists. In the last case, one may want to find a partial solution optimizing certain criteria, for example, as many satisfied constraints as possible. We restrict our discussion in finding one solution. B. Definition of constraint satisfaction agents An agent used to solve CSPs is represented as a constraint satisfaction agent (CSAgent), and is defined as follows: Definition 1: A constraint satisfaction agent, a, represents an element in the search space, S, with energy being equal to (3) ∀a∈S, Energy(a) = − ¦im=1 χ (a, Ci )

1 a violates Ci . The purpose of each where χ (a, Ci ) =  0 otherwise CSAgent is to maximize the energy by the behaviors it can take.

SMCB-E-03172004-0141.R2

3

As can be seen, the energy of each CSAgent is the negative value of the number of violated constraints. Therefore, the higher the energy is, the better the CSAgent is. When the energy of a CSAgent increases to 0, the CSAgent becomes a solution. The domain for each variable is finite and discrete, thus the elements can be numbered by natural numbers, that is to say, the domain D={d1, d2, …, d|D|} is equivalent to D={1, 2, …, |D|}. When all domains are transformed to the sets of natural numbers, the solutions of some CSPs take on specific characteristics. On the basis of them, CSPs can be divided into two types: Definition 2: Let S be the set of solutions for a CSP and all domain sets be represented by the sets of natural numbers. If ∀s∈S is a permutation of 1, 2, …, n, the CSP is defined as a permutation CSP; otherwise it is defined as a non-permutation CSP. For a non-permutation CSP, the search space is still S. Since permutation CSPs are a special case of non-permutation CSPs, the search space can be reduced to the set of all permutations of 1, 2, …, n, that is, S P = { P1 , P2 , ..., Pn !} and Pi = 〈 Pi1 , Pi 2 , ..., Pi n 〉, 1 ≤ i ≤ n ! and Pi j ∈ {1, 2, ..., n} , 1 ≤ j ≤ n

(4) When one uses EAs to solve problems, the search space must be encoded such that individuals can be represented in a uniform encoding. For example, GAs usually encode the search space by the binary coding, thus, each individual is a bit sequence. For permutation CSPs, it is natural to represent each individual as a permutation of 1, 2, …, n. The search space is SP with size of n!. For non-permutation CSPs, the most straightforward encoding method is to represent each individual as an element of S, which is used by most existing algorithms. Under this method, the search space is |D1|×|D2|×…×|Dn|. Meanwhile, some methods use a permutation coding with a corresponding decoder. For example, in [27], [28], each individual is represented by a permutation of the variables, and the permutation is transformed to a partial instantiation by a simple decoder that considers the variables in the order they occur in the permutation and assigns the first possible domain value to that variable. If no value is possible without introducing a constraint violation, the variable is left uninstantiated. In what follows, this decoder is labeled as Decoder1 and all permutations of the variables is labeled as SPx . In fact, P1

)

, xP2 , , xPn ⌡ ∈ SPx  ( ∑ P1 , P2 , , Pn ⌡ ∈ S P )

S Px = {∑ x1 , x2 , x3 ⌡, ∑ x1 , x3 , x2 ⌡, ∑ x2 , x1 , x3 ⌡,

∑ x2 , x3 , x1 ⌡, ∑ x3 , x1 , x2 ⌡, ∑ x3 , x2 , x1 ⌡}

.

(6)

According to Decoder1, each element in SPx can be transformed to a partial instantiation of the variables, namely,  x1 , x2 , x3 → 1, 1, * x1 , x3 , x2 → 1, 1, *  (7) x2 , x3 , x1 → 1, *, 1  x2 , x1 , x3 → 1, 1, *  x3 , x2 , x1 → 1, *, 1  x3 , x1 , x2 → 1, 1, * where the “*” represents that the corresponding variable is left uninstantiated. As can be seen, no element in SPx can be  transformed to a solution for the CSP by Decoder1. For CSAgents, we want to design an encoding method that not only covers solutions for any CSPs, but also maintains a manageable search space. Therefore, on the basis of Decoder1, we propose minimum conflict encoding (MCE). In MCE, each CSAgent is represented as not only an element of SPx , but also an element of S, so that some behaviors can be designed to deal with the values of the variables directly. As a result, each CSAgent must record some information, and is represented by the following structure:

 Pk ≠ Pl and ∀k , l ∈ {1, 2, ..., n} , ( k ≠ l ) ⇒  k l  Pi ≠ Pi

( ∀∑ x

Example 2: The CSP is given in Example 1 with the following SPx

(5)

Therefore, SPx is equivalent to SP, and the size of the search space is also n!. Decoder1 uses a greedy algorithm, thus a serious problem exists. That is, for some CSPs, Decoder1 cannot decode any permutation to a solution. As a result, the algorithms based on Decoder1 may not find solutions at all. Here is a simple example:

CSAgent = Record P: A permutation, P∈SP for permutation CSPs, and P ∈ SPx for non-permutation CSPs; V: V∈S, the value for each variable, and for non-permutation CSPs only; E: The energy of the CSAgent, E=Energy(P) for permutation CSPs, and E=Energy(V) for non-permutation CSPs; SL: The flag for the self-learning behavior, which will be defined later. If SL is True, the self-learning behavior can be performed on the CSAgent; otherwise cannot; End. In the following text, CSAgent(•) is used to represent the corresponding component in the above structure. By MCE, each CSAgent includes both a permutation and an assignment to all variables. In fact, the assignment is what we need. So minimum conflict decoding (MCD) is designed corresponding to MCE, which transforms P to V. The main idea of MCD is to assign the value violating minimum number of constraints to a variable. The variables are considered in the order they occur in the permutation, and only the constraints between the variable and those previously assigned values are taken into account. After MCD is performed, no variable is left uninstantiated. The details of MCD are shown in Algorithm 1.

Algorithm 1

Minimum conflict decoding

SMCB-E-03172004-0141.R2

4

Input: CSAgent: the CSAgent needs to be decoded; Pos: the position to start decoding; Output: CSAgent(V); Let CSAgent ( P ) = ∑ xP1 , xP2 , , xPn ⌡ , CSAgent(V)=〈v1, v2, …, vn〉, and Conflicts (vi ) = ∑ j =1 χ (vi , C j ) m

(8)

1 vi violates C j where χ (vi , C j ) =  . Conflicts(vi) only 0 otherwise considers the variables assigned values. MinC and MinV are the current minimum number of conflicts and the corresponding value. begin if (Pos=1) then begin vP1 := 1 ;

i := 2; end else i := Pos; repeat vPi := 1 ;

MinC := Conflicts(vPi ) ; MinV := 1; j := 2; repeat vPi := j ; if ( Conflicts(vPi ) < MinC ) then begin MinC := Conflicts(vPi ) ;

MinV := j; end; j := j+1; until ( j >| DPi | );

vPi := MinV ; i := i+1; until (i>n); end. Two kinds of behaviors are designed for CSAgents, namely, the behaviors performing on P (P-behavior) and the behaviors performing on V (V-behavior). When V-behaviors are performed on a CSAgent, the energy can be updated directly. But when P-behaviors are performed on a CSAgent, it must be decoded by MCD before updating the energy. If MCD starts to decode from the first variable in the permutation for any CSAgent, the information generated by V-behaviors would loss. Therefore, we set a parameter, Pos, for MCD, which is the first position changed by P-behaviors. Thus, the value of the variables before Pos is left untouched such that some

information generated by V-behaviors can be preserved. For a new CSAgent, Pos is set to 1. MCD(CSAgent, Pos) represents MCD is performed on CSAgent. In fact, Algorithm 1 is just a general implementation of the idea of MCD, and can be used to all kinds of CSPs. However, for specific CSPs, the idea of MCD can be combined with the knowledge of the problems, with a more efficient implementation obtained. For example, the degree of each vertex can be used when dealing with GCPs, and the details will be shown in Section IV.B.1. C. Environment of constraint satisfaction agents In order to realize the local perceptivity of agents, the environment is organized as a latticelike structure, which is similar to our previous work in [34]. Definition 3: All CSAgents live in a latticelike environment, L, which is called an agent lattice. The size of L is Lsize×Lsize, where Lsize is an integer. Each CSAgent is fixed on a lattice-point and can only interact with the neighbors. Suppose that the CSAgent located at (i, j) is represented as Li,j, i, j=1,2,…,Lsize, then the neighbors of Li,j, Neighborsi,j, are defined as follows: (9) Neighborsi , j = { Li ′, j , Li , j ′ , Li′′, j , Li , j ′′ } i + 1 i ≠ Lsize i − 1 i ≠ 1  j −1 j ≠ 1 , , j′ =  , i ′′ =  where i ′ =  1 1 L i = L j = i = Lsize  size  size 1

 j + 1 j ≠ Lsize . j ′′ =  j = Lsize 1 The agent lattice can be represented as the one in Fig.1. Each circle represents a CSAgent, the data represent the position in the lattice, and two CSAgents can interact with each other if and only if there is a line connecting them. In traditional EAs, those individuals used to generate offspring are usually selected from all individuals according to their fitness. Therefore, the global fitness distribution of a population must be determined. But in nature, a global selection does not exist, and the global fitness distribution cannot be determined either. In fact, the real natural selection only occurs in a local environment, and each individual can only interact with those around it. That is, in each phase, the natural evolution is just a kind of local phenomenon. The information can be shared globally only after a process of diffusion. In the aforementioned agent lattice, to achieve their purposes, CSAgents will compete with others so that they can gain more resources. Since each CSAgent can only sense the local environment, the behaviors can only take place between the CSAgent and the neighbors. There is no global selection at all, so the global fitness distribution is not required. A CSAgent interacts with the neighbors so that information is transferred to them. In such a manner, the information is diffused to the whole agent lattice. As can be seen, the model of the agent lattice is closer to the real evolutionary mechanism in nature than the model of the population in traditional EAs.

SMCB-E-03172004-0141.R2 D. Behaviors of constraint satisfaction agents The purpose of an algorithm for solving CSPs is to find solutions at a low computational cost as much as possible. So the computational cost can be considered as the resources of the environment in which all CSAgents live. Since the resources are limited and the behaviors of the CSAgents are driven by their purposes, a CSAgent will compete with others to gain more resources. On the bases of this, three behaviors are designed for CSAgents to realize their purposes, that is, the competitive behavior, the self-learning behavior, and the mutation behavior. The former two belong to P-behaviors, whereas the last belongs to V-behaviors. Competitive behavior: In this behavior, the energy of a CSAgent is compared with those of the neighbors. The CSAgent can survive if the energy is maximum; otherwise the CSAgent must die, and the child of the one with maximum energy among the neighbors will take up the lattice-point. Suppose that the competitive behavior is performed on the CSAgent located at (i, j), Li,j, and Maxi,j is the CSAgent with maximum energy among the neighbors of Li,j, namely, Maxi,j∈Neighborsi,j and ∀CSAgent∈Neighborsi,j, then CSAgent(E)≤Maxi,j(E). If Li,j(E) ≤ Maxi,j(E), then Maxi,j generates a child CSAgent, Childi,j, to replace Li,j, and the method is shown in Algorithm 2; otherwise Li,j is left untouched.

Algorithm 2 Competitive behavior Input: Maxi,j: Maxi,j(P)=〈m1, m2, …, mn〉 for permutation CSPs, and Maxi , j ( P ) = ∑ xm1 , xm2 , , xmn ⌡ for non-permutation CSPs; pc: A predefined parameter, and a real number between 0-1; Output: Childi,j: Childi,j(P)=〈c1, c2, …, cn〉 for permutation CSPs, and Childi , j ( P ) = ∑ xc1 , xc2 , , xcn ⌡ for non-permutation CSPs; Swap(x, y) exchanges the values of x and y. U(0, 1) is a uniform random number between 0 and 1. Random(n, i) is a random integer among 1, 2, …, n and is not equal to i. Min(i, j) is the smaller one between i and j. begin Childi,j(P) := Maxi,j(P); i := 1; Pos := n+1; repeat if (U(0, 1) 0 (27) ij

tBest

Therefore, ∀k≥i, pi.k≥pij.k>0.  It follows from this theorem that there is always a positive probability to transit from an agent lattice to the one with identical or higher energy and a zero probability to the one with lower energy. Thus, once MAEA-CSPs enters L0, it will never escape. Theorem 5: Multiagent evolutionary algorithm for constraint satisfaction problems converges to the global optimum as time tends to infinity. Proof: It is clear that one can consider each Li, i∈E, as a state in a homogeneous finite Markov chain. According to theorem 4, the transition matrix of the Markov chain can be written as follows: 0  0   p0.0   p p  0  C 0  -1.-1 = P ′ =  -1.0 (28)        R T    p- m.0 p- m.-1  p- m.-m  Obviously, R=(p-1.0 p-2.0 … p-m.0)T>0, T≠0, C=(p0.0)=(1)≠0. According to theorem 2, P′∞ is given by,  Ck 0   ∞ 0 ∞ k  = C (29) P ′ = lim P ′ = lim k -1 i  ∞  k −i k  k →∞ k →∞  R 0 ∑ T RC T     i =0  where C∞=(1), R∞=(1 1 … 1)T. Thus, P′∞ is a stable stochastic

SMCB-E-03172004-0141.R2 matrix, and 1 0  0    1 0  0 P ′∞ =        1 0  0    Therefore, lim Pr{Energy ( Lt ) = 0} = 1 t →∞

9

(30)

(31)

where Pr stands for the probability. This implies that multiagent evolutionary algorithm for constraint satisfaction problems converges to the global optimum as time tends to infinity.  IV. EXPERIMENTAL STUDIES ON NON-PERMUTATION CONSTRAINT SATISFACTION PROBLEMS Reference [54] indicates that any CSP can be equivalently transformed to a binary CSP. Since binary CSPs belong to non-permutation CSPs, they are used to test the performance of MAEA-CSPs in the section. Systematic analyses are also made to check the effect of the parameters, pc and pm. Afterwards, MAEA-CSPs is used to solve a more practical case, graph coloring problems. Four parameters must be assigned before MAEA-CSPs is used to solve problems. Lsize×Lsize is equivalent to the population size in traditional EAs, so Lsize is generally selected from 3 to 10, and is set to 5 in the following experiments. pc and pm control the competitive behavior and the mutation behavior, and are set to 0.2 and 0.05. EvaluationMax is set to 105 for binary CSPs, and 104 for GCPs. All experiments are executed on a 2.4-GHz Pentium IV PC with 1G RAM. A. Binary constraint satisfaction problems The experiments are conducted on the test suite [78] used to compare the performances of 11 available algorithms in [29]. The suite consists of 250 solvable problem instances. Each instance has 20 variables, and the domain set is D={Di | Di={1, 2, …, 20} and 1≤i≤20}. The 250 instances are divided into 10 groups according to their difficulty, p={0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33}, and each group has 25 instances. The higher the value of p, the more difficult the instance is. We perform 10 independent runs on each of the 25 instances belonging to a given p value, amounting to 250 data points used for calculating our measures. The method used to measure the performance of MAEA-CSPs is similar to that of [29], namely, the effectiveness and the efficiency are considered. The effectiveness is measured by the success rate (SR) and the mean error (ME). The success rate is the percentage of runs finding a solution. Since all instances in the test suite are solvable, the highest value of the SR is 100%. The error is defined for a single run as the number of constraints not satisfied by the best CSAgent in the agent lattice when the run terminates. For any given set of runs, the ME is the average of these error values. This measure provides information on the quality of partial solutions. This can be a useful way of comparing algorithms having equal SR values lower than 100%. The efficiency is measured by the average number of

evaluations to solution (AES), which has been widely used to measure the computational cost of EAs. For a given algorithm, it is defined as the average number of evaluations of successful runs. Consequently, when the SR=0, the AES is undefined. Note also, when the number of successful runs is relatively low, the AES is not statistically reliable. Therefore, it is important to consider the SR and AES together to make a clear interpretation of the results. Among the three measures, the most important is the success rate since we ultimately want to solve problems. The second is the AES, in the case of comparable SR figures, we prefer the faster EA. The ME is only used as a third measure as it might suggest preferences between EAs failed. A.1 Comparison between MAEA-CSPs and the existing algorithms Reference [29] has made a comparison among 11 available algorithms, and the results show that the best four algorithms are H-GA.1 [20], [21], H-GA.3 [20], [21], SAW [27], [28], and Glass-Box [19]. Therefore, we reimplemented the four algorithms according to the descriptions and parameters given in [29], and made a comparison between them and MAEA-CSPs. The results of the four algorithms are consistent with those of [29]. The comparison results are shown in Fig.2. In order to facilitate other researchers to make a comparison between MAEA-CSPs and their algorithms, the results of MAEA-CSPs are also shown in Table I. As can be seen, the SR of MAEA-CSPs is the highest among the five algorithms with 100% SR for p=0.24-0.26. Since the ME of SAW is much larger than those of the other algorithms, the ME of SAW and MAEA-CSPs is compared in top-right of Fig.2(b). Fig.2(b) shows that the ME of MAEA-CSPs is also the best among the five algorithms. Since the AES is not statistically reliable when the SR is very low, only the AES of the instances whose SR is larger than 10% is plotted in Fig.2(c). All the SR of the five algorithms for p=0.30-0.33 is very low. But the ME of MAEA-CSPs is small, and only 1 to 3 constraints are not satisfied. To summarize, MAEA-CSPs outperforms the four other algorithms. A.2 Parameter analyses of MAEA-CSPs In this experiment, pc and pm are increased from 0.05 to 1 in steps of 0.05, with 20×20=400 groups of parameters obtained. MAEA-CSPs with the 400 groups of parameters is used to solve the 10 groups of instances. According to the SR, the 10 groups of instances can be divided into 3 classes. The first class includes the instances with low p values, that is, p=0.24-0.26. Since the SR for this class is higher than 90% for most of the parameters and the ME is very low, the graphs for the SR and the AES are shown in Fig.3. The second class includes the instances with p varying from 0.27 to 0.29. Since the SR for this class is lower, the graphs for the SR and the ME are shown in Fig.4. The last class includes the instances with high p values. Since the instances in this class are very difficult, only the graphs for the ME are shown in Fig.5. As can be seen, pc has a larger effect on the performance of MAEA-CSPs. For the first class, the SR is higher than 90%

SMCB-E-03172004-0141.R2

10

when pc is larger than 0.1. Although the AES increases with pc when pc is larger than 0.2, the AES is smaller when pc is in 0.1-0.3. For the second and the last class, the results are similar, namely, when pc is in 0.1-0.3, the SR is higher, and the ME is smaller. Although pm does not affect the performance of MAEA-CSPs obviously, Fig.4 and Fig.5 show the SR is slightly higher and the ME is slightly smaller when pm is small. To summarize, it is better to choose pc from 0.1-0.3 and pm from 0.05-0.3. In addition, although the performance of MAEA-CSPs with above parameters is better, MAEA-CSPs still performs stably when pc is larger than 0.2. It shows that the performance of MAEA-CSPs is not sensitive to the parameters, and MAEA-CSPs is quite robust and easy to use. B. Graph coloring problems Many problems of practical interest can be modeled as graph coloring problems. The purpose is to color each vertex of a graph by a certain color, such that no two vertices incident to any edge have the same color. Given an undirected graph G={V, E}, where V={V1, V2, …, Vn} is the set of vertices, and E={ei,j | Vi and Vj are connected} is the set of edges. An m-coloring problem can be described as follows when representing each vertex as a variable, and the domain corresponds to the given set of m colors:   x = {x1 , x2 ,  , xn }, xi represents Vi ∈ V   D = {D1 , D2 ,  , Dn },  Di = {1, 2,  , m}, i = 1, 2,  , n (32)   C = C ({xi , x j }) = 〈d , d 〉 | ∀i, j ∈ {1, 2, ..., n},  ∀d ∈ {1, 2, ..., m}, ei , j ∈ E}  In what follows, MCD for GCPs is described first, and then MAEA-CSPs is tested on the official benchmarks from DIMACS graph coloring challenge [79]. Finally, a comparison is made between MAEA-CSPs and two recent algorithms.

{

B.1 Minimum conflict decoding for graph coloring problems Since the constraints only occur between two connected vertices, and the degree of a vertex is far smaller than |V| in general, we can record the degree and the list of vertices connected to an identical vertex in advance, and only check the lists when performing MCD. Algorithm 5

Minimum conflict decoding for graph coloring problems Input: CSAgent: the CSAgent needs to be decoded; Pos: the position to start decoding; Output: CSAgent(V); Let CSAgent ( P ) = ∑ xP1 , xP2 ,  , xPn ⌡ , CSAgent(V)=〈v1,

v2, …, vn〉, degreePi vertex

xPi ,

stands for the degree of

and

connected to xPi .

vertex Pji

the

jth

vertex

begin if (Pos=1) then begin vP1 := 1 ;

i := 2; end else i := Pos; repeat for j:=1 to m do Conflictsj:=0; for j:=1 to degreePi do if ( vertex Pji then begin

has been assigned value)

Color := the color of vertex Pji ; ConflictsColor := ConflictsColor + 1; end; Color := j, where Conflictsj is minimum among Conflicts1, Conflicts2, …, Conflictsm; vPi := Color ; i := i+1; until (i>n); end. In Algorithm 5, the second “for” only considers the vertices connected to xPi and previously assigned values. When this “for” scans once, the number of vertices with each color are computed out. So the algorithm is more efficient. B.2 Experimental results of MAEA-CSPs Presently, the DIMACS challenge has 79 problems. These problems include various graphs, such as random graphs (DSJCx.y), geometric random graphs (DSJRx.y), “quasirandom” graphs (flatx_y_0), graphs from real-life applications (fpsol2.i.x, inithxi.x, mulsol.i.1, and zeroin.i.1). To investigate the performance of MAEA-CSPs extensively, all the 79 problems are considered, and the results averaged over 10 independent runs are shown in Table II. In the DIMACS challenge, some problems counted each edge once whereas some ones counted twice. Since the graphs are undirected, |E| in Table II counts each edge once. The value of m is set to the number of optimal colors when known, and the smallest value we can find in the literature (values in parentheses) when unknown. In Table II, “c” stands for the number of constraints not satisfied. Since each constraint corresponds to an edge, “c” is the number of edges whose two vertices have the same color. As can be seen, MAEA-CSPs finds exact solutions in all 10 runs for 33 out of 79 problems, and c of these problems is 0. Moreover, the running times of these problems are very short. “(1-c/|E|)×100” indicates the percentage of edges colored properly. Among the 79 problems, there are 66 problems that MAEA-CSPs colors more than 99% of the edges properly, 9 problems that 98-99% of the edges properly, 2 problems that 97-98% of the edges properly, and 2 problems that 94-95% of

SMCB-E-03172004-0141.R2 the edges properly. The running times of MAEA-CSPs are shorter than 1 second for 46 problems, between 1 to 4 seconds for 22 problems, between 5 to 10 seconds for 3 problems, and longer than 10 seconds for only 8 problems. B.3 Comparison between MAEA-CSPs and two recent algorithms Reference [55] proposed an algorithm, labeled as SSC, for GCPs based on energy functions and neural networks. Reference [55] tested SSC by a part of the DIMACS challenge, and good performance was obtained. Table III compares SSC with MAEA-CSPs, where the results of SSC are those reported in [55]. SSC counted each violated edge twice, which has been confirmed by the author of [55]. So we correct the “c” of SSC in Table III. MAEA-CSPs finds better solutions than SSC for 74 out of 79 problems. Reference [33] has also designed an agent evolutionary model, ERA, and ERA was tested on 11 problems out of the DIMACS challenge. Table IV compares ERA with MAEA-CSPs, where the results of ERA are those reported in [33]. What should be noted is [33] gives the number of vertices having no conflicts, which is labeled as “v” in Table IV, whereas MAEA-CSPs gives the number of edges having conflicts. For all the 11 problems, both the running times of ERA and of MAEA-CSPs are 0.00s. Moreover, MAEA-CSPs colors 100% of the edges properly for these 11 problems, that is to say, 100% of the vertices are colored properly. So the performance of MAEA-CSPs is better than that of ERA. In conclusion, Table III and Table IV show that MAEA-CSPs outperforms the two algorithms. V. EXPERIMENTAL STUDIES ON PERMUTATION CONSTRAINT SATISFACTION PROBLEMS A. n-queen problems n-queen problems are typical problems of CSPs, and have been studied in depth [20], [33], [56]-[59]. A solution to an n-queen problem places n queens on an n×n chessboard such that no two queens can attack each other. The purpose of this study is twofold. First, a general search algorithm to n-queen problems is a useful indicator of the potential for solving other CSPs [57]. Second, an n-queen problem is itself a model of the maximum coverage problems [58]. A solution to n-queen problems guarantees that each object can be accessed by the outside world from any one of the eight neighboring directions (two vertical, two horizontal, and four diagonal directions) without a conflict with other objects. This result can directly apply to VLSI testing, air traffic control, modern communication systems, data/message routing in a multiprocessor computer. In what follows, the description of n-queen problems and an example for the self-learning behavior are given first, and then the experimental results. A.1 Problem description In any solution of n-queen problems, no two queens are on

11 the same row or on the same column, because it would cause collisions. Therefore, a solution for n-queen problems can be represented by a permutation of 1, 2, …, n, P=〈P1, P2, …, Pn〉. It denotes to put the ith queen on the ith row, the Pith column. Thus, n-queen problems belong to permutation CSPs, and the search space is SP. When we encode each candidate solution by a permutation, the collisions on rows and columns are avoided naturally, and what left is to find the permutations satisfying the constraints on the diagonal lines. An n×n grid has (2n-1) positive diagonal lines and (2n-1) negative diagonal lines, and they have the following characteristics: the difference between the row index and the column index is constant on any positive diagonal lines, and the sum of both indexes is constant on any negative diagonal lines (see Fig.6). Consequently, an n-queen problem can be described as follows:  x = {x1 , x2 ,  , xn }   D = {D1 , D2 ,  , Dn },  Di = {1, 2,  , n}, i = 1, 2, , n  (33)  C = C ({xi , x j }) = 〈d k , d l 〉 | ∀i, j , d k , dl ∈ {1, 2,  , n},  (d k = d l ) or (i − j = d k − d l ) or   (i − j = d l − d k )}

{

To explain the self-learning behavior explicitly, an example of the performing process is given. Example 3: Let the problem under consideration be an 8-queen problem, and the CSAgent on which the self-learning behavior performed be 〈1, 2, 8, 4, 5, 7, 3, 6〉. Then the performing process is given in Fig.7. For more clarity, Fig.7 is explained further. Since the queen in the first row has collisions, the algorithm first searches for a swap, see Fig.7(b). Suppose that the selected l is 3, then after the swap is performed, we have Energynew=-6>Energyold=-7. So the swap is successful, Repeat is set to True, and Iteration increases 1. Here the queen in the first row has no collisions, so the algorithm goes on to deal with the queen in the second row, see Fig.7(c). Suppose that 3 is chosen for l. Although the swap is successful, the queen in the second row still has collisions. Therefore, the algorithm still deals with it, see Fig.7(d). Suppose that 4 is chosen for l, then the swap is successful and the algorithm goes on to deal with the queen in the third row, see Fig.7(e). Suppose that 5,7,1,2 are chosen for l in turn, but all swaps failed. At the moment, Iteration is equal to 7, so Iteration is set to 1 and the algorithm goes on to deal with the queen in the fourth row, see Fig.7(f). Suppose that 6 is chosen for l, then the swap is successful. Presently, the queen in the fourth row has no collisions, thus the algorithm goes on to deal with the queens in the following rows, see Fig.7(g). Since the queens in the sixth and eighth row have no collisions, the algorithm deals with the queens in the fifth and seventh rows. But all swaps failed. However, Repeat is True, hence the algorithm restarts from the beginning. During this time, the collisions of the queens in the fifth and seventh rows cannot be eliminated. At which point Repeat is equal to False and the algorithm halts. The final state of the CSAgent is shown in Fig.7(h) and the energy increases

SMCB-E-03172004-0141.R2 from –7 to –1.

12



A.2 Experimental results for n-queen problems There exist solutions for n-queen problems with n greater than or equal to 4 [58]. So the termination criterion of MAEA-CSPs is set to find a solution, namely, set EvaluationMax to ∞. Thus, all the following experimental results have SR=100% and ME=0. Then, what should be studied further is the computational cost, which is measured by the running times. Reference [20] designed a GA that can solve 104-queen problems, but the result was reported in a figure, and cannot be compared with that of MAEA-CSPs directly. The agent evolutionary model, ERA [33], can solve 7000-queen problems. Therefore, the results of MAEA-CSPs for 103-104-queen problems are given first, and then a comparison is made between MAEA-CSPs and ERA [33]. The results are shown in Table V. The results of ERA are obtained by running the software [80] on the same environment. The software restrains n in 4-7000. As can be seen, the performance of MAEA-CSPs is much better than that of ERA, and MAEA-CSPs only uses 9 milliseconds to solve the 104-queen problems. In order to study the scalability of MAEA-CSPs along the problem scale, MAEA-CSPs is used to solve 5×104-107-queen problems. In this experiment, n increases from 5×104 to 107 in steps of 5×104. At each sampled value of n, 50 trials are carried out, and the average running times are shown in Fig.8. As can be seen, the running times of MAEA-CSPs can be approximated by the function, (5.04×10-6×n1.07). That is to say, MAEA-CSPs has a linear time complexity for n-queen problems. For more clarity, the average running times and the standard deviation of MAEA-CSPs are shown in Table VI for the problems with 1×106, 2×106, …, 1×107 queens. MAEA-CSPs only uses 13s to solve 106-queen problems, and 150 seconds to solve 107-queen problem. Moreover, all standard deviations are very small, and the maximum is only 1.15. So MAEA-CSPs not only has a fast convergence rate, but also has a stable performance. B. Job-shop scheduling problems Since modern manufacturing environments are very complex, making it very difficult and time-consuming for people to create good schedules, it is a great advantage to have the scheduling process performed automatically by a computer system. There has been a considerable research effort on scheduling, such as the methods based on EAs [60]-[65], multiagents [66], simulated annealing [67], neural networks [68], hybrid heuristic technique [69], fuzzy logic [70]-[72]. The focus of this section is on job-shop scheduling problems [73]. A small modification has been made on MAEA-CSPs so that JSPs can be solved, which is introduced in Section V.B.1. Then, the experimental results are reported in Section V.B.2. B.1 Multiagent evolutionary algorithm for job-shop scheduling problems A JSP of size n×m consists of n jobs and m machines. For each job Ji, a sequence of m operations Oi=(oi,1, oi,2, …, oi,m)

describing the processing order of the operations of Ji is given. Each operation oi,j is to be processed on a specific machine and has a processing time τi,j. When the operations are processed each machine can process only one operation at a time, each job can only have one operation processed at a time, and no preemption can take place. A solution to a JSP is a schedule specifying when to process each of the operations, not violating any of the constraints. One encoding method in common use for JSPs is permutation with repetition, where a schedule is described as a sequence of all n×m operations, and each operation in the sequence is described by the job-number. CSAgent(P) is represented as this format, namely, P = P1 , P2 , ..., Pn×m and ( P (1) = m) and (34) ( P (2) = m) and ... and ( P (n) = m) Where Pi∈{1, 2, …, m}, i=1, 2, …, n×m, and P(j), j=1, 2, …, n stand for the number of j in P. When transforming P into a schedule, Pi stands for oj,k if Pi=j and the number of j among P1, P2, …, Pi is equal to k. The schedule is obtained by considering the operations in the order they occur in P and assigning the earliest allowable time to that operation. The energy of a CSAgent is equal to the negative value of the makespan, where the makespan of a schedule is defined as the time elapsed from the beginning of processing until the last job has finished. Such encoding method has the advantage that no infeasible solutions can be represented, and each CSAgent corresponds to a feasible schedule. Therefore, the problem is transformed to a constraint optimization problem in fact. Although MAEA-CSPs is designed for constraint satisfaction problems, it can also deal with JSPs by making a small modification. Since a number in P occurs many times, the algorithm must be prevented from swapping two identical values. So the sixth line in Algorithm 2, “l := Random(n, i)”, should be changed to “l := Random(n, i) and ci≠cl”.

Let P=〈P1, P2, …, Pn×m〉,

o pi represent the operation

corresponding to Pi, and M pi be the machine on which o pi is to be processed. Suppose that Pi and Pj (Pi≠Pj and i