**4. An algorithm for the construction of covering arrays using a simulated annealing technique**

Often the solution space of an optimization problem has many local minima. A simple local search algorithm proceeds by choosing random initial solution and generating a neighbor from that solution. The neighboring solution is accepted if it is a cost decreasing transition. Such a simple algorithm has the drawback of often converging to a local minimum. The simulated annealing algorithm (SA), though by itself it is a local search algorithm, avoids getting trapped in a local minimum by also accepting cost increasing neighbors with some probability. SA is a general-purpose stochastic optimization method that has proven to be an effective tool for approximating globally optimal solutions to many types of NP-hard combinatorial optimization problems. In this section, we briefly review SA algorithm and propose an implementation to solve CAC problem.

SA is a randomized local search method based on the simulation of annealing of metal. The acceptance probability of a trial solution is given by Eq. 3, where *T* is the *temperature* of the system, Δ*E* is the difference of the costs between the trial and the current solutions (the cost change due to the perturbation), Eq. 3 means that the trial solution is accepted by nonzero probability *e*(−Δ*E*/*T*) even though the solution deteriorates (*uphill move*).

$$\mathcal{E}(P) = \begin{cases} 1 & \text{if } \Delta E < 0 \\ \mathcal{e}^{( - \frac{\Delta E}{T})} & \text{otherwise} \end{cases} \tag{3}$$

Uphill moves enable the system to escape from the local minima; without them, the system would be trapped into a local minimum. Too high of a probability for the occurrence of uphill moves, however, prevents the system from converging. In SA, the probability is controlled by temperature in such a manner that at the beginning of the procedure the temperature is sufficiently high, in which a high probability is available, and as the calculation proceeds the temperature is gradually decreased, lowering the probability (Jun & Mizuta, 2005).

#### **4.1 Internal representation**

The following paragraphs will describe each of the components of the implementation of our SA. The description is done given the matrix representation of an CA. An CA can be represented as a matrix *M* of size *N* × *k*, where the columns are the parameters and the rows are the cases of the test set that is constructed. Each cell *mi*,*<sup>j</sup>* in the array accepts values from the set {1, 2, ..., *vj*} where *vj* is the cardinality of the alphabet of *j th* column.

#### **4.2 Initial solution**

The *initial solution M* is constructed by generating *M* as a matrix with maximum Hamming distance. The Hamming distance *d*(*x*, *y*) between two rows *x*, *y* ∈ *M* is the number of elements in which they differ. Let *ri* be a row of the matrix *M*. To generate a random matrix *M* of maximum Hamming distance the following steps are performed:

1. Generate the first row *r*<sup>1</sup> at random.

10 Grid Computing

parallelizing our SA algorithm. Section 6 describes how to implement our algorithm on a

Often the solution space of an optimization problem has many local minima. A simple local search algorithm proceeds by choosing random initial solution and generating a neighbor from that solution. The neighboring solution is accepted if it is a cost decreasing transition. Such a simple algorithm has the drawback of often converging to a local minimum. The simulated annealing algorithm (SA), though by itself it is a local search algorithm, avoids getting trapped in a local minimum by also accepting cost increasing neighbors with some probability. SA is a general-purpose stochastic optimization method that has proven to be an effective tool for approximating globally optimal solutions to many types of NP-hard combinatorial optimization problems. In this section, we briefly review SA algorithm and

SA is a randomized local search method based on the simulation of annealing of metal. The acceptance probability of a trial solution is given by Eq. 3, where *T* is the *temperature* of the system, Δ*E* is the difference of the costs between the trial and the current solutions (the cost change due to the perturbation), Eq. 3 means that the trial solution is accepted by nonzero

<sup>1</sup> *i f*Δ*<sup>E</sup> <* <sup>0</sup>

Uphill moves enable the system to escape from the local minima; without them, the system would be trapped into a local minimum. Too high of a probability for the occurrence of uphill moves, however, prevents the system from converging. In SA, the probability is controlled by temperature in such a manner that at the beginning of the procedure the temperature is sufficiently high, in which a high probability is available, and as the calculation proceeds the

The following paragraphs will describe each of the components of the implementation of our SA. The description is done given the matrix representation of an CA. An CA can be represented as a matrix *M* of size *N* × *k*, where the columns are the parameters and the rows are the cases of the test set that is constructed. Each cell *mi*,*<sup>j</sup>* in the array accepts values from

The *initial solution M* is constructed by generating *M* as a matrix with maximum Hamming distance. The Hamming distance *d*(*x*, *y*) between two rows *x*, *y* ∈ *M* is the number of elements in which they differ. Let *ri* be a row of the matrix *M*. To generate a random matrix *M* of

*<sup>T</sup>* ) *otherwise* (3)

*th* column.

probability *e*(−Δ*E*/*T*) even though the solution deteriorates (*uphill move*).

*e*(<sup>−</sup> <sup>Δ</sup>*<sup>E</sup>*

temperature is gradually decreased, lowering the probability (Jun & Mizuta, 2005).

(*P*) =

the set {1, 2, ..., *vj*} where *vj* is the cardinality of the alphabet of *j*

maximum Hamming distance the following steps are performed:

**4. An algorithm for the construction of covering arrays using a simulated**

grid architecture.

**annealing technique**

**4.1 Internal representation**

**4.2 Initial solution**

propose an implementation to solve CAC problem.


$$\log(r\_l) = \sum\_{s=1}^{i-1} \sum\_{v=1}^{k} d(m\_{\text{s,p}}, m\_{i,\text{\textnu}}), \text{where } d(m\_{\text{s,p}}, m\_{i,\text{\textnu}}) = \begin{cases} 1 \text{ if } m\_{\text{s,p}} \neq m\_{i,\text{\textnu}} \\ 0 \text{ otherwise} \end{cases} \tag{4}$$

4. Repeat from step 2 until *M* is completed.

An example is shown in Fig. 2; the number of symbols different between rows *r*<sup>1</sup> and *c*<sup>1</sup> are 4 and between *r*<sup>2</sup> and *c*<sup>1</sup> are 3 summing up 7. Then, the hamming distance for the candidate row *c*<sup>1</sup> is 7.

$$\begin{array}{c} \text{Rows} \begin{cases} r\_1 = \begin{Bmatrix} 2 \ 1 \ 0 \ 1 \end{Bmatrix} \\ r\_2 = \begin{Bmatrix} 1 \ 2 \ 0 \ 1 \end{Bmatrix} \\ c\_1 = \begin{Bmatrix} 0 \ 2 \ 1 \ 0 \end{Bmatrix} \end{array} \end{array} \qquad \begin{cases} d(r\_1, c\_1) = 4 \\ d(r\_2, c\_1) = 3 \\ g(c\_1) = 7 \end{Bmatrix}$$

Fig. 2. Example of the hamming distance between two rows *r*1, *r*<sup>2</sup> that are already in the matrix M and a candidate row *c*1.

#### **4.3 Evaluations function**

The *evaluation function E*(*M*) is used to estimate the goodness of a candidate solution. Previously reported metaheuristic algorithms for constructing CA have commonly evaluated the quality of a potential solution (covering array) as the number of combination of symbols missing in the matrix *M* (Cohen et al., 2003; Nurmela, 2004; Shiba et al., 2004). Then, the expected solution will be zero missing.

In the proposed SA implementation this evaluation function definition was used. Its computational complexity is equivalent to *O*(*N*( *k t* )).

#### **4.4 Neighborhood function**

Given that our SA implementation is based on Local Search (LS) then a neighborhood function must be defined. The main objective of the neighborhood function is to identify the set of potential solutions which can be reached from the current solution in a LS algorithm. In case two or more neighborhoods present complementary characteristics, it is then possible and interesting to create more powerful compound neighborhoods. The advantage of such an approach is well documented in (Cavique et al., 1999). Following this idea, and based on the results of our preliminary experimentations, a neighborhood structure composed by two different functions is proposed for this SA algorithm implementation.

Two *neighborhood functions* were implemented to guide the local search of our SA algorithm. The neighborhood function N1(*s*) makes a random search of a missing *t*-tuple, then tries by setting the *j th* combination of symbols in every row of *<sup>M</sup>*. The neighborhood function <sup>N</sup>2(*s*) randomly chooses a position (*i*, *j*) of the matrix *M* and makes all possible changes of symbol. During the search process a combination of both N1(*s*) and N2(*s*) neighborhood functions is employed by our SA algorithm. The former is applied with probability *P*, while the latter is employed at an (1 − *P*) rate. This combined neighborhood function N3(*s*, *x*) is defined in Eq. 5, where *x* is a random number in the interval [0, 1).

$$\mathcal{N}\_3(s, \mathbf{x}) = \begin{cases} \mathcal{N}\_1(s) & \text{if } \mathbf{x} \le P \\ \mathcal{N}\_2(s) & \text{if } \mathbf{x} > P \end{cases} \tag{5}$$

**Algorithm 4:** Sequential simulated annealing for the CAC problem

*<sup>T</sup>* ) *> x* **then**

**<sup>3</sup> repeat**

**<sup>4</sup> for** *i* ← 1 **to** *L* **do**

**<sup>12</sup> end if <sup>13</sup> end if <sup>14</sup> end for**

described as follows:

**<sup>3</sup> repeat**

**<sup>7</sup> end if**

participating processes.

**<sup>5</sup> if** *E*(*M*) *< E*(*M�*) **then**

**<sup>8</sup>** *CALCULATE*\_*CONTROL*(*T*,*φ*) **<sup>9</sup> until** *termination condition is satisfied*;

back to step 2, else algorithm terminates.

**Algorithm 5:** Simulated annealing for the master node

**<sup>8</sup> if** <sup>Δ</sup>*<sup>E</sup> <* <sup>0</sup> **or** *<sup>e</sup>*(<sup>−</sup> <sup>Δ</sup>*<sup>E</sup>*

**<sup>10</sup> if** *E*(*M*) *< E*(*M�*) **then**

**<sup>15</sup>** *CALCULATE*\_*CONTROL*(*T*,*φ*) **<sup>16</sup> until** *termination condition is satisfied*;

**<sup>1</sup>** *INITIALIZE*(*M*, *T*, *L*) ; /\* Create the *initial solution*. \*/ **<sup>2</sup>** *<sup>M</sup>�* <sup>←</sup> *<sup>M</sup>* ; /\* Memorize the *best solution*. \*/

Using Grid Computing for Constructing Ternary Covering Arrays 233

**<sup>5</sup>** *Mi* ← *GENERATE*(*M*); /\* Perturb current state. \*/ **<sup>6</sup>** Δ*E* ← *E*(*Mi*) − *E*(*M*) ; /\* Evaluate cost function. \*/ **<sup>7</sup>** *x* ← *random* ; /\* Range [0,1). \*/

**<sup>9</sup>** *M* ← *Mi* ; /\* Accept new state. \*/

**<sup>11</sup>** *<sup>M</sup>�* <sup>←</sup> *<sup>M</sup>* ; /\* Memorize the *best solution*. \*/

process will probably get a different solution at the end of iterations. This approach may be

1. The *master* node set *T* = *T*0, generates an *initial*\_*solution* using the Hamming distance

4. If the termination condition is not met, each process reduces the temperature and goes

Algorithm 5 shows the pseudocode for *master* node. The function INITIALIZE computes a start solution (using Hamming distances algorithm) and initial values of the parameters *T* and *L*. The *master* node distributes the initial parameters to slave nodes, and awaits the results. Each *L* iterations, the slaves send their results to the master node (See Algorithm 6). The master node selects the best solution. If the termination criterion is not satisfied, the master node computes a new value for the parameter *T* (cooling schedule) and the number of

**<sup>1</sup>** *INITIALIZE*(*M*, *T*, *L*) ; /\* Create the *initial solution*. \*/ **<sup>2</sup>** *<sup>M</sup>�* <sup>←</sup> *<sup>M</sup>* ; /\* Memorize the *best solution*. \*/

**<sup>4</sup>** *<sup>M</sup>* <sup>←</sup> *annealing*\_*worker*(*M�*, *<sup>T</sup>*, *<sup>L</sup>*) ; /\* Call workers \*/

**<sup>6</sup>** *<sup>M</sup>�* <sup>←</sup> *<sup>M</sup>* ; /\* Memorize the *best solution*. \*/

2. At the current temperature *T*, each worker begins to execute iterative operations (*L*). 3. At the end of iterations, the *master* is responsible for collecting the solution obtained by each process at current temperature and broadcasts the best solution of them among all

algorithm (See Section 4.2) and distributes them to each workers.

consecutive temperature decrements with no improvement in the solution.

In the next section it is presented our parallel simulated annealing approach for solving CAC problem.

### **4.5 Cooling schedule**

The *cooling schedule* determines the degree of uphill movement permitted during the search and is thus critical to the SA algorithm's performance. The parameters that define a cooling schedule are: an initial temperature, a final temperature or a stopping criterion, the maximum number of neighboring solutions that can be generated at each temperature, and a rule for decrementing the temperature. The literature offers a number of different cooling schedules, see for instance (Aarts & Van Laarhoven, 1985; Atiqullah, 2004). In our SA implementation we preferred a geometrical cooling scheme mainly for its simplicity. It starts at an initial temperature *Ti* which is decremented at each round by a factor *α* using the relation *Tk* = *<sup>α</sup>Tk*−<sup>1</sup> . For each temperature, the maximum number of visited neighboring solutions is *<sup>L</sup>*. It depends directly on the parameters (*N*, *k* and *v* is the maximum cardinality of *M*) of the studied covering array. This is because more moves are required for CAs with alphabets of greater cardinality.

#### **4.6 Termination condition**

The *stop criterion* for our SA is either when the current temperature reaches *Tf* , when it ceases to make progress, or when a valid covering array is found. In the proposed implementation a lack of progress exists if after *φ* (frozen factor) consecutive temperature decrements the best-so-far solution is not improved.

#### **4.7 SA Pseudocode**

The Algorithm 4 presents the simulated annealing heuristic as described above. The meaning of the four functions is obvious: INITIALIZE computes a start solution and initial values of the parameters *T* and *L*; GENERATE selects a solution from the neighborhood of the current solution, using the neighborhood function N3(*s*, *x*); CALCULATE\_CONTROL computes a new value for the parameter *T* (cooling schedule) and the number of consecutive temperature decrements with no improvement in the solution.

#### **5. Parallel simulated annealing**

In this section we propose a parallel strategy to construct CAs using a simulated annealing algorithm.

A common approach to parallelizing simulated annealing is to generate several perturbations in the current solution simultaneously. Some of them, those with small variance, locally explore the region around the current point, while those with larger variances globally explore the feasible region. If each process has got different perturbation or move generation, each 12 Grid Computing

is employed at an (1 − *P*) rate. This combined neighborhood function N3(*s*, *x*) is defined in

In the next section it is presented our parallel simulated annealing approach for solving CAC

The *cooling schedule* determines the degree of uphill movement permitted during the search and is thus critical to the SA algorithm's performance. The parameters that define a cooling schedule are: an initial temperature, a final temperature or a stopping criterion, the maximum number of neighboring solutions that can be generated at each temperature, and a rule for decrementing the temperature. The literature offers a number of different cooling schedules, see for instance (Aarts & Van Laarhoven, 1985; Atiqullah, 2004). In our SA implementation we preferred a geometrical cooling scheme mainly for its simplicity. It starts at an initial temperature *Ti* which is decremented at each round by a factor *α* using the relation *Tk* = *<sup>α</sup>Tk*−<sup>1</sup> . For each temperature, the maximum number of visited neighboring solutions is *<sup>L</sup>*. It depends directly on the parameters (*N*, *k* and *v* is the maximum cardinality of *M*) of the studied covering array. This is because more moves are required for CAs with alphabets of

The *stop criterion* for our SA is either when the current temperature reaches *Tf* , when it ceases to make progress, or when a valid covering array is found. In the proposed implementation a lack of progress exists if after *φ* (frozen factor) consecutive temperature decrements the

The Algorithm 4 presents the simulated annealing heuristic as described above. The meaning of the four functions is obvious: INITIALIZE computes a start solution and initial values of the parameters *T* and *L*; GENERATE selects a solution from the neighborhood of the current solution, using the neighborhood function N3(*s*, *x*); CALCULATE\_CONTROL computes a new value for the parameter *T* (cooling schedule) and the number of consecutive temperature

In this section we propose a parallel strategy to construct CAs using a simulated annealing

A common approach to parallelizing simulated annealing is to generate several perturbations in the current solution simultaneously. Some of them, those with small variance, locally explore the region around the current point, while those with larger variances globally explore the feasible region. If each process has got different perturbation or move generation, each

<sup>N</sup>1(*s*) if *<sup>x</sup>* <sup>≤</sup> *<sup>P</sup>*

<sup>N</sup>2(*s*) if *<sup>x</sup> <sup>&</sup>gt; <sup>P</sup>* (5)

Eq. 5, where *x* is a random number in the interval [0, 1).

problem.

**4.5 Cooling schedule**

greater cardinality.

**4.7 SA Pseudocode**

algorithm.

**4.6 Termination condition**

best-so-far solution is not improved.

**5. Parallel simulated annealing**

decrements with no improvement in the solution.

N3(*s*, *x*) =


**Algorithm 4:** Sequential simulated annealing for the CAC problem

process will probably get a different solution at the end of iterations. This approach may be described as follows:


Algorithm 5 shows the pseudocode for *master* node. The function INITIALIZE computes a start solution (using Hamming distances algorithm) and initial values of the parameters *T* and *L*. The *master* node distributes the initial parameters to slave nodes, and awaits the results. Each *L* iterations, the slaves send their results to the master node (See Algorithm 6). The master node selects the best solution. If the termination criterion is not satisfied, the master node computes a new value for the parameter *T* (cooling schedule) and the number of consecutive temperature decrements with no improvement in the solution.

**Algorithm 5:** Simulated annealing for the master node


(sequential) or more (parallel) computing nodes. In the case that no free computing nodes are available, jobs are queued. Thus, the load of a CE must be considered when estimating

Using Grid Computing for Constructing Ternary Covering Arrays 235

3. WN (Working Nodes): Each one of the computing resources accessible through a CE. Due to the heterogeneous nature of Grid infrastructure, the response time of a job will depend

4. SE (Storage Element): Storage resources in which a task can store long-living data to be used by the computers of the Grid. This practice is necessary due to the size limitation imposed by current Grid Middlewares in the job file attachment (10 Megabytes in the gLite case). So, use cases which require the access to files which exceed that limitation are forced to use these Storage Elements. Nevertheless, the construction of CAs is not a data-intensive

5. LFC (Logic File Catalog): A hierarchical directory of logical names referencing a set of

6. BDII (Berkley Database Information System): Service point for the Information System which registers, through LDAP, the status of the Grid. Useful information relative to CEs,

7. R-GMA (Relational Grid Monitoring Architecture). Service for the registration and

8. VOMS (Virtual Organisation Management System). Authorization infrastructure to define

A production infrastructure such as EGI involves tens of thousands of resources from hundreds of sites, involving tens of countries and a large human team. Since it is a general-purpose platform, and although there is a common middleware and a recommended operating system, the heterogeneity in the configuration and operation of the resources is inevitable. This heterogeneity, along with other social and human factors such as the large geographical coverage and the different skills of operators introduces a significant degree of uncertainty in the infrastructure. Even considering that the service level required is around 95%, it is statistically likely to find in each large execution sites that are not working properly. Thus, prior to beginning the experiments, it is necessary to do empirical tests to define a group of valid computing resources (CEs) and this way facing resource setup problems. These tests can give some real information like computational speed, primary and secondary memory sizes and I/O transfer speed. These data, in case there are huge quantities of resources, will

Once the computing elements, where the jobs will be submitted, have been selected, the next step involves correctly specifying the jobs. In that sense, it will be necessary to produce the specification using the job description language in gLite. An example of a JDL file can be seen

the turnaround of a job.

on the characteristics of the WN hosting it.

use case and thus the use of SEs can be avoided.

WNs and SEs can be obtained by querying this element.

notification of information in most of the EGI services.

**6.2 Preprocessing task: Selecting the most appropriate CEs**

be helpful to establish quality criteria choosing resources.

Some of this terms will be referenced along the following sections.

physical files and replicas stored in the SEs.

the access rights to resources.

**6.3 Asynchronous schema**

in the Fig. 3.

#### **Algorithm 6:** Worker algorithm

