Danny Munera

Daniel Diaz

Salvador Abreu

Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization

Keywords: QAP, extremal optimization, heuristics, parallelism, cooperation

Several real-life applications can be stated in terms of the Quadratic Assignment Problem. Finding an optimal assignment is computationally very difficult, for many useful instances. We address this problem using a local search technique, based on Extremal Optimization and present experimental evidence that this approach is competitive. Moreover, cooperative parallel versions of our solver improve performance so much that large and hard instances can be solved quickly.

Introduction

The Quadratic Assignment Problem (QAP) was introduced in 1957 by Koopmans and Beckmann [START_REF] Koopmans | Assignment Problems and the Location of Economic Activities[END_REF] as a model of a facilities location problem. This problem consists in assigning a set of n facilities to a set of n specific locations minimizing the cost associated with the flows of items among facilities and the distance between them. This combinatorial optimization problem has many other real-life applications: scheduling, electronic chipset layout and wiring, process communications, turbine runner balancing, data center network topology, to cite but a few [START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF]. Unfortunately this problem is known to be NP-hard and finding efficient algorithms to solve it has attracted a lot of research in recent years.

Exact (or complete) methods like dynamic programming, cutting plane techniques and branch & bound procedures have been successfully applied to mediumsize QAP instances but cannot solve larger instances (e.g. when n > 20). To tackle these problems, one must resort to incomplete methods which are designed to quickly provide good, albeit sub-optimal, solutions. This is the case of approximation algorithms, i.e. algorithms running in polynomial time yet able to guarantee solutions within a constant factor of the optimum. Unfortunately, it is known that there is no -approximation algorithm for the QAP [START_REF] Sahni | P-Complete Approximation Problems[END_REF]. Another class of incomplete methods is provided by meta-heuristics. Since 1990 several meta-heuristics have been successfully applied to the QAP: tabu search, simulated annealing, genetic algorithms, GRASP, ant-colonies [START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF]. The current trend is to specialize existing heuristics, to compose different meta-heuristics (hybrid procedures) and to use parallelism.

In this paper we propose EO-QAP: an Extremal Optimization (EO) procedure for QAP. EO is a nature-inspired general-purpose meta-heuristics to solve combinatorial optimization problems [START_REF] Boettcher | Nature's way of optimizing[END_REF]. This local search procedure, a priori, has several advantages: it is easy to implement, it does not get confounded by local minima and takes only one adjustable parameter. We experimentally demonstrate that EO-QAP performs well on the set of QAPLIB benchmark instances. It is, however, known that it is difficult with EO to have fine control on the trade-off between search intensification and diversification: some strategies have been proposed to overcome this limitation [START_REF] Randall | Intensification Strategies for Extremal Optimisation[END_REF], but they entail a more complex tuning process. In this paper we put forth two other approaches which contribute to a more effective handling of QAP using EO: firstly, we propose a simple extension to the original EO which allows the user to have more control over the stochastic behavior of the algorithm. Secondly, we propose to use cooperative parallelism to promote more intensification and/or diversification. Our implementation uses a parallel framework [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF] written in X10 [START_REF] Charles | X10: An Object-Oriented Approach to Non-Uniform Cluster Computing[END_REF][START_REF] Saraswat | A Brief Introduction to X10 (for the High Performance Programmer)[END_REF]. We show that the cooperative parallel version behaves very well on the hardest instances.

The rest of the paper is organized as follows. Section 2 discusses QAP and provides the necessary background. Section 3 presents our EO algorithm for QAP and proposes an extension to the original EO. Several experimental results are laid out and discussed in section 4 and we conclude in Section 5.

Background

Before introducing the main object of this paper, we need to recall some background topics: the Quadratic Assignment Problem (QAP), Extremal Optimization (EO) and Cooperative Parallel Local Search (CPLS).

QAP

Since its introduction in 1957, QAP has been widely studied and several surveys are available [START_REF] Burkard | Quadratic Assignment Problems[END_REF][START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Loiola | A survey for the quadratic assignment problem[END_REF][START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF].

A QAP problem of size n consists of two n × n matrices (a ij ) and (b ij ). Let Π(n) be the set of all permutations of {1, 2, . . . n}, the goal of QAP is to find a permutation π ∈ Π(n) which minimizes the following objective function:

F (π) = n i=1 n j=1 a ij • b πiπj (1)
For instance, in facility location problems, the a matrix represents interfacility flows and b encodes the inter-location distances. In that context, both matrices are generally symmetric: ∀ i,j a ij = a ji and b ij = b ji . However, in other settings the matrices can become asymmetric. Indeed, QAP can be used to model scheduling, chip placement and wiring on a circuit board, to design typewriter keyboards, for process communications, for turbine runner balancing among many other applications [START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Zaied | A Survey of Quadratic Assignment Problems[END_REF]. The computational difficulty of QAP stems form the fact that the objective function contains products of variables (hence the term quadratic) and in the fact that the theoretical search space of an instance of size n is the set of all permutations Π(n) whose cardinality is n!. In 1976, Sahni and Gonzalez proved that QAP is NP-hard [START_REF] Sahni | P-Complete Approximation Problems[END_REF] (the famous traveling salesman problem can be formulated as a QAP). Moreover, the same authors proved that there is no -approximation algorithm for QAP (unless P=NP). In practice QAP is one of the most difficult combinatorial optimization problems with many real-life applications.

QAP can be (optimally) solved with exact methods like dynamic programming, cutting plane techniques and branch & bound algorithms (together with efficient lower bound methods). Constraint Programming does not work well on QAP and, surprisingly, SAT solvers have not been extensively used for QAP. However, general problems of medium size (e.g. n > 20) are out of reach for these methods (even if some particular larger instances can be solved). It is thus natural to use heuristics to solve QAP. In the last decades several meta-heuristics were successfully applied to QAP: tabu search, simulated annealing, genetic algorithms, GRASP, ant-colonies [START_REF] Said | A Comparative Study of Meta-heuristic Algorithms for Solving Quadratic Assignment Problem[END_REF]. In this paper we propose to use Extremal Optimization to attack QAP problems.

Extremal Optimization

In 1999, Boettcher and Percus proposed the Extremal Optimization (EO) procedure [START_REF] Boettcher | Nature's way of optimizing[END_REF][START_REF] Boettcher | Extremal Optimization: an Evolutionary Local-Search Algorithm[END_REF][START_REF] Boettcher | Extremal Optimization[END_REF] as a meta-heuristics to solve combinatorial optimization problems. EO is inspired by self-organizing processes often found in nature. It based on the concept of Self-Organized Criticality (SOC) initially proposed by Bak [START_REF] Bak | Self-Organized Crtiticality: An Explenation of 1/f Noise[END_REF][START_REF] Bak | How Nature Works: The Science of Self-organized Criticality[END_REF] and in particular by the Bak-Sneppen model of SOC [START_REF] Bak | Punctuated equilibrium and criticality in a simple model of evolution[END_REF]. In this model of biological evolution, species have a fitness ∈ [0, 1] (0 representing the worst degree of adaptation). At each iteration, the species with the worst fitness value is updated, i.e. its fitness is replaced by a new random value. This change also affects all other species connected to this "culprit" element and their fitness value also gets updated. This results in an extremal process which progressively eliminates the least fit species (or forces them to mutate). Repeating this process eventually leads to a state where all species have a good fitness value, i.e. a SOC. EO follows this line: it inspects the current configuration (assignment of variables), selects the worst variable (the one having the lowest fitness) and replaces its value by a random value. However, always selecting the worst variable can lead to a deterministic behavior and the algorithm can stay blocked in a local minimum. To avoid this, the authors propose an extended algorithm which first ranks the variables in increasing order of fitness (the worst variable has thus a rank k = 1) and then resorts to a Probability Distribution Function (PDF) over the ranks k in order to introduce uncertainty in the search process:

P (k) = k -τ (1 ≤ k ≤ n) (2)
This power-law probability distribution takes a single parameter τ which is problem-dependent. Depending on the value of τ , EO provides a wide variety of search strategies from pure random walk (τ = 0) to deterministic (greedy) search (τ → ∞). With an adequate value for τ , EO cannot be trapped in local minima since any variable is susceptible to mutate (even if the worst are privileged). This parameter can be tuned by the user. Moreover, the original paper proposes a default value depending on n: τ = 1 + 1 ln(n) . EO displays several a priori advantages: it is a simple meta-heuristic (it can be easily programmed), it is controlled by only one free parameter (a fine tuning of several parameters becomes quickly tedious) and it does not need to be aware about local minima. Nevertheless, EO has been successfully applied to large-scale optimization problems like graph bi-partitioning, graph coloring, Spin Glasses or the traveling salesman problem [START_REF] Boettcher | Extremal Optimization: an Evolutionary Local-Search Algorithm[END_REF]. Boettcher and Percus point out, however, that depending on the problem, "a drawback of the EO method is that a general definition of fitness for the individual variables may prove ambiguous or even impossible" [START_REF] Boettcher | Optimization with Extremal Dynamics[END_REF]. To overcome this, De Sousa and Ramos proposed an extension called Generalized Extremal Optimization [START_REF] De Sousa | Function optimization using extremal dynamics[END_REF][START_REF] De Sousa | Generalized Extremal Optimization for Solving Complex Optimal Design Problems[END_REF]. Zhou and al. proposed a variant called Continuous Extremal Optimization to deal with continuous optimization problems [START_REF] Zhou | Continuous Extremal Optimization for Lennard-Jones Clusters[END_REF]. It has been also argued that one main issue with EO is that is does not provide a fine control of the intensification. Randall and Lewis propose some intensification strategies to improve EO [START_REF] Randall | Intensification Strategies for Extremal Optimisation[END_REF]. We present two alternatives to overcome this limitation: we propose a simple extension to improve the stochastic capabilities of EO and we show how cooperative parallelism can help to achieve intensification or diversification through communications.

Cooperative Parallel Local Search

Parallel local search methods have been proposed in the past [START_REF] Alba | Parallel Metaheuristics: A New Class of Algorithms[END_REF][START_REF] Alba | Parallel Metaheuristics: Recent Advances and New Trends[END_REF][START_REF] Diaz | Parallel Constraint-Based Local Search on the Cell/BE Multicore Architecture[END_REF]. In this article we are interested in multi-walk methods (also called multi-start) which consist in a concurrent exploration of the search space, either independently or cooperatively via communication between processes. The Independent Multi-Walks method (IW) [START_REF] Verhoeven | Parallel Local Search[END_REF] is the easiest to implement since the solver instances do not communicate with each other. However, the resulting gain tends to flatten when scaling over a hundred of processors [START_REF] Caniou | Large-scale parallelism for constraint-based local search: the costas array case study[END_REF], and can be improved upon. In the Cooperative Multi-Walks (CW) method [START_REF] Toulouse | Systemic Behavior of Cooperative Search Algorithms[END_REF], the solver instances exchange information (through communication), hoping to hasten the search process. However, implementing an efficient cooperative method is a very complex task: several choices have to be made about the communication which influence each other and which are problem-dependent [START_REF] Toulouse | Systemic Behavior of Cooperative Search Algorithms[END_REF].

We build on the framework for Cooperative Parallel Local Search (CPLS) proposed in [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF][START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. This framework, available as an open source library in the X10 programming language, allows the programmer to tune the search process through an extensive set of parameters which, in the present version, statically condition the execution. CPLS augments the IW strategy with a tunable communication mechanism, which allows for the cooperation between the multiple instances to seek either an intensification or diversification strategy for the search. At present, the tuning process is done manually: we have not yet experimented with parameter self-adaptation (still an experimental feature).

The basic component of CPLS is the explorer node which consists in a local search solver instance. The point is to use all the available processing units by mapping each explorer node to a physical core. Explorer nodes are grouped into teams, of size N P T (see Figure 1). This parameter is directly related to the trade-off between intensification and diversification. N P T can take values from 1 to the maximum number of nodes (frequently linked to maximum number of available cores in the execution). When N P T is equal to 1, the framework coincides with the IW strategy, it is expected that each 1-node team be working on a different region of the search space, without seek parallel intensification. When N P T is equal to the maximum number of nodes (creating only 1 team in the execution), the framework has the maximum level of parallel intensification, but it is not possible to enable parallel diversification between teams. Each team seeks to intensify the search in the most promising neighborhood found by any of its members. The parameters which guide the intensification are the Report Interval (R) and Update Interval (U ): every R iterations, each explorer node sends its current configuration and the associated cost to its head node. The head node is the team member which collects and processes this information, retaining the best configurations in an Elite Pool (EP ) whose size |EP | is parametric. Every U iterations, explorer nodes randomly retrieve a configuration from the EP , in the head node. An explorer node may adopt the configuration from the EP , if it is "better" than its own current configuration with a probability pAdopt. Simultaneously, the teams implement a mechanism to cooperatively diversify the search, i.e. they try to extend the search to different regions of the search space.

Typically, each problem benefits from intensification and diversification on some level. Therefore, the tuning process of the CPLS parameters seeks to provide the appropriate balance between the use of the intensification and diversification mechanisms, in hope of reaching better performance than the noncooperative parallel solvers (e.g. Independent Multi-Walks). A detailed description of this framework may be found in [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF].

3 EO-QAP: an EO Procedure for QAP

General Procedure

Our EO-QAP algorithm starts from a random permutation π ∈ Π(n), with an associated cost given by F (π). To ensure we only consider proper permutations, we only perform swap operations on pairs of elements from any given one: this is how value assignment is classically implemented in permutation problems, thereby eschewing the costly explicitly encoding of an all-different constraint. We define the permutation resulting from swapping π i and π j :

π i↔j = µ | µ i = π j , µ j = π i , µ k = π k ∀k ∈ {i, j} (3)
The neighborhood of π i is the set N (π i ) of the permutations obtained from π by swapping π i with any another value. By extension, the neighborhood of π is the set N (π) of all permutations obtained by swapping any two values:

N (π i ) = {π i↔j | 1 ≤ j ≤ n, j = i} (4)
N (π) = n i=1 N (π i ) (5)
Most local search procedures (in particular hill climbing) would select, among all elements of N (π), the best neighbor µ, i.e. the one minimizing the next global cost function. By doing so, they have to deal with the problem of being trapped in a local minimum. Instead, EO defines a fitness value λ i for each value π i , with the understanding that a value with a low fitness is more likely to mutate, i.e. to get swapped. We define the fitness value λ i as the best possible improvement of the cost F when moving to a π i 's neighbor.

λ i = min µ∈N (πi) F (µ) -F (π) (6)
A negative λ thus represents an improvement of the cost. At each iteration, the n fitness values are evaluated and ranked, with rank k = 1 for the worst fitness. EO-QAP will thus favor the mutation of a value which improves (decreases) the objective function. The value to mutate is chosen stochastically from a probability distribution over the rank order. This comes down to pick a value at rank k (1 ≤ k ≤ n) with a probability P (k) = k -τ . Let π r be the value measured by the k th fitness value, then π r will be deemed the "culprit" and be forced to mutate. For this we need to choose a target value π s for the swap. Several possibilities exist: one is to choose π s randomly. Another, as in the original EO article, is to pick a random value using the same probability distribution. We propose a third possibility applying the min-conflict heuristic [START_REF] Minton | Minimizing Conflicts: A Heuristic Repair Method for Constraint-Satisfaction and Scheduling Problems[END_REF]: select the best possible value, that is, the value which minimizes the objective function of the next configuration. The algorithm then swaps π r and π s (thus moving to a neighbor) and iterates with this new configuration. The process stops when a some condition is reached (e.g. a time limit or a given cost is reached). The best solution found so far is then returned (see Algorithm 1).

Algorithm 1 EO-QAP: an EO procedure for QAP

Implementation Commentary and Complexity Analysis.

A permutation π can be encoded by an array of n integers sol[] (with sol[i] = π i ). Line 6 computes the fitness λ i for each value π i . Thanks to Taillard and his famous Robust Taboo Search (RoTS), we know how to compute the evolution of the global cost after a swap incrementally, instead of recomputing it each time from scratch -see equations (1) and (2) which define ∆(µ, ., .) in [START_REF] Taillard | Comparison of iterative searches for the quadratic assignment problem[END_REF]. This results in an evaluation of all λ in O(n 2 ), while a naïve algorithm is in O(n 3 ). The simplest data structure to manage λ is just an array fit[] whose elements are pairs of the form <index,lambda>. Initially we have fit[i]=<i, λ i >. Line 7 sorts the fit[] array on the second field (λ), in ascending order, this can be done in O(n log 2 (n)). Line 8 picks a value at position k with a probability P (k). Since the PDF and τ are constant along the execution of the algorithm, it is more efficient to pre-compute the n samples of P (k), (1 ≤ k ≤ n) and store then in an array prob[]. To pick a random k, we can use a roulette-wheel selection on prob[] in O(n) theoretically (but closer to O(1) in practice due to the PDF). It is also possible to use a a binary search in O(log 2 (n)) storing the cumulative PDF (prob[k] = k i=1 P (i)). Line 9: the variable to mutate is given by sol[fit [k].index]. Line 10 selects the other variable for the swap, as per the min-conflict heuristic; this is done in O(n). However, as this variable had already been found when computing λ (Line 6), one just has to record it. To this end, the elements of the fit[] array are refactored as <index,lambda,index2>, where index2 contains the index of the variable which minimizes the cost. The second variable for the swap is then simply given by sol[fit [k].index2]. The overall complexity of each iteration (i.e. of the main loop body) is thus O(n 2 ).

Extending Extremal Optimization

Because it relies on just one parameter, EO is comparatively very simple to use (tuning many local search parameters can become very laborious). As we show in section 4, the results of EO-QAP are very good on many instances of QAPLIB. However, some harder instances need a fine control of the tradeoff between intensification and diversification. EO handles these two strategies with the same tool: the probability distribution P (k) = k -τ . Depending on the returned value, either the current path continues to be improved (with a high probability) or, in the extreme case, completely abandoned (with a low probability). Every variable has a non-zero chance of being selected and EO is not affected by local extrema. The choice of the probability distribution P has thus a great impact, also determined by its parameter τ , which must be selected by the user. Tuning τ for some hard problems (e.g. tai40a) turned out to be difficult. We therefore decided to extend EO so as to accept different probability distribution functions (PDF). The user can thus choose the most appropriate PDF. For simplicity, all proposed PDFs take a single input parameter τ , the other parameters (if any) being either constant or functionally dependent on τ . We now discuss a few interesting PDFs, and how they are influenced by the single τ parameter:

PDF definition

Usage in EO Name The EO algorithm behaves very differently, depending on the PDF and the chosen value for τ . Figure 2 shows the curves associated with these PDFs for a size n = 40, picking different values for τ . Clearly, with the power law, the first ranked variable has a very high probability to be selected, the probability then decreases very fast but the variables with a high rank (i.e."good" variables) are very likely to mutate (e.g. when τ = 0.5). This results in less intensification, and suits some problems perfectly. If this behavior is unwanted, a larger value for τ may be used, nevertheless, this rapidly puts too strong a pressure on the besk ranked variables which, for some problems, will be selected too frequently. It can be difficult to find a good trade-off. On the other hand, the exponential law skews probabilities a little bit more in favor of lower ranked variables. This is clear from the shapes of the PDFs (see Figure 2). The normal (Gaussian) law is also interesting because the curve decreases slowly for the very first ranks, then more rapidly and then more slowly again: the first ranked "worst" variables will then be selected with a high but comparable priority. Finally, we found the gamma law interesting because it is not strictly decreasing and can, for instance, give more priority to the second or third variable than to the first one. With this PDF we obtained good results for tai40a with τ = 1.5. It is worth noticing that a shifted normal law can also be used as a non-stricly decreasing function, e.g. using P (k) = N ormal(k, 2, e τ ). It is worth noticing that the best PDF and τ combination is not the same when EO-QAP is run sequentially or in parallel.

P ower(x, τ ) = x -τ P (k) = P ower(k, τ ) power law Expon(x, µ) = µ e -µx P (k) = Expon(k, τ ) exponential law N ormal(x, µ, σ) = 1 σ √ 2π e -(x-µ) 2 2σ 2 P (k) = N ormal(k, 1, τ ) normal law Gamma(x, k, θ) = 1 Γ (k)θ k x k-1 e -x θ P (k) = Gamma(k,
Clearly, allowing different PDFs enhances the power of EO which can attack more problems efficiently. The user now has more precise control over the behavior of the algorithm, which remains simple with only two tunable parameters: the PDF and its τ value. It is even possible allow the user to provide his own customized PDF as a file of n values P 1 , . . . P n . At run-time, the prob[] array above mentionned is populated as follows: prob[k] = P k n i=1 Pi (each value being divided by the sum of all values to ensure the whole PDF = 1).

Experimental Evaluation

In this section we present experimental results on the entire QAPLIB test set. To do this, we developed an X10 implementation of EO-QAP. 3 Because EO is a stochastic procedure, we ran each problem 10 times and averaged the results. The number of possible experiments is very high: with different PDFs and τ values, varying the timeout, testing sequential or parallel runs, with different topologies and communication strategies, etc. We thus adopted a 3-stage protocol:

1. Attempt to solve all QAPLIB problems (134 instances) with a basic version of EO-QAP (i.e. without any tuning) and a very short timeout (in order to be able to try 10 runs for each of the 134 instances). All problems for which the Best Known Solution (BKS) was reached for every execution are definitively classified as solved. The others (solved less than 10 times) form the test set of the next stage.

2. To attack the remaining problem we ran EO-QAP (with same parameters and timeout) in independent parallel multi-walks, without communication, on 32 cores of a single machine. As previously, we collected the fully solved instances which need no longer be considered. The remaining instances form the input set for the next stage. 3. The remaining problems are the hardest ones: for these, we used a cooperative parallel version of EO-QAP on 128 cores, tuning the PDF and τ value, and using a larger timeout of 10 minutes.

Stage 1: Sequential Execution

In this first stage, the input test set consists of the 134 QAPLIB instances, which are run sequentially on an AMD Opteron 6376 clocked at 2.3 GHz, using a single core. This is the basic version of EO-QAP with the original power-law PDF and default value for τ (see Section 2.2). For each problem we report the BKS (which is sometimes the optimum), the number of times the BKS is reached (#BKS), i.e. the number of times the problem is solved, the average execution time (in seconds) and the the Average Percentage Deviation (APD) which is the average of the 10 relative deviation percentages computed as follows: 100 F (sol)-BKS

BKS

. We use a short timeout of 5 minutes. Even if the solver stops as soon as the BKS is reached a limited timeout is needed to be able to run the 134 instances 10 times.

Table 1 presents the whole results. Surprisingly, even with this straightforward and suboptimal setting, EO-QAP performs quite well: more than 50% of the instances get totally solved. More precisely: 68 instances are fully solved (in green), and among the remaining 66, only 25 are never solved. On average, the 41 others are solved 4.6 times (in orange). The average APD for the 66 instances not fully solved is about 2.2%.

Stage 2: Independent Parallelism

In this stage, we ran the EO-QAP algorithm in parallel without communication with the same settings as in the first stage (default PDF, default τ , timeout 5 min). The machine was the same: a quad AMD Opteron 6376 clocked at 2.3 GHz, but using 32 cores. These parameters make it possible to assess what improvement we can easily obtain by means of parallel execution. Indeed, this from of parallelism (sometimes called embarrassingly parallelism) works by performing multiple independent walks to explore the search space. Each worker blindly explores a region of the search space, looking for a solution. The process ends as soon as any solver reaches a solution. Since all EO solvers start from a random point, we can expect that they will all visit different regions of the search space (i.e. ensuring a form of diversification), thus increasing the chance to find a solution. Such a parallelization of an algorithm is easy to implement and often behaves very well (see Section 2.3).

The results of this experiment are summarized in Table 2. This form of parallelism brings a significant improvement in performance and reach. Exactly already provides library functions for all these operations. The resulting solver is then composed of several EO-QAP instances running in parallel, which cooperate by communicating in order to converge faster on a solution. As per [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF], the CPLS parameters which control the cooperation are as follows:

-Team Size (N P T ): we tested various configurations and defined N P T = 16.

There are thus 8 teams composed of 16 explorer nodes each running the EO-QAP procedure. This is constant for all problems. -Report Interval (R): we manually tuned it (starting from the average number of iterations collected during the previous stage divided by 10). -Update Interval (U ): we experimented different ratio and retained U/R = 2.

-Elite Pool (EP ): its size is fixed to 4 for all problems.

-pAdopt: is set to 1. An EO-QAP instance receiving a better configuration than its current one always switches to this new one.

This experiment has been carried out on a cluster of 16 machines, each with 4 × 16-core AMD Opteron 6376 CPUs running at 2.3 GHz and 128 GB of RAM. The nodes are interconnected with InfiniBand FDR 4× (i.e. 56 GBPS.) We had access to 4 nodes and used up to 32 cores per node, i.e. 128 cores. We stay with a timeout of 5 minutes. Finally, we tested deeply two PDFs (power and exponential) and their τ value and retained the best combination for each problem instance (we could not yet test the Normal law and the Gamma law, while efficient in sequential seems, not well suited for parallelism or else using a different τ value).

Table 3 presents the results obtained on the hardest instances. The table also reports the average number of iterations, the report and update interval (R and U ), the number of times the winning algorithm has adopted an elite configuration, the PDF and τ value used. In this last stage, 15 new problems become fully solved. Moreover the average time to solve them is only 24.5 seconds. Only 8 remain unsolved. On average, the remaining 10 get solved 5 out of 10 times. Moreover, the average APD over the 18 not fully solved is only 0.25%. Even when the optimum is not reached the returned solution is close to this optimum.

The table shows that an efficient execution corresponds to a limited number of "adoptions" of an elite configuration (less than 5 changes). This is directly correlated to the value of R and U . There are some exceptions like tai25b and tho40 which are both fully solved. We plan to analyze in details these both situations (e.g. varying R and U).

Regarding the power and exponential PDF, there is no winner. It is worth noticing that sometimes the difference is huge. For instance, sko90 is solved 9 times with the power law but only 2 times with the exponential law. The reverse occurs for taiXX a for which the exponential law performs much better.

The performance of the cooperative parallel EO-QAP, is on par with the best competing solutions, while retaining a much simpler internal structure [START_REF] James | A Cooperative Parallel Tabu Search Algorithm for the Quadratic Assignment Problem[END_REF]. Due to space limitations, we do not develop this further.

Conclusion and Further Work

We have proposed EO-QAP: an Extremal Optimization (EO) procedure for QAP. The basic sequential version of EO-QAP, while simple behaves rather well on several instances of QAPLIB. To attack hardest instances we first proposed a simple extension to the original EO procedure allowing for different probability distribution functions (PDF). The user can select the most adequate PDF depending on the degree of intensification wanted. Moreover, we resorted to cooperative parallelism using a framework written in the X10 parallel language, which provides the user with a fine degree of control of the intensification and diversification for the search. The cooperative version of EO-QAP displays very good results: using 128 cores, 116 instances of QAPLIB are systematically solved at each execution, 10 instances are solved half the time and only 8 instances remain unsolved (but the obtained solution is near to the optimum).

We now plan to attack other known hard instances. Future work includes the study of a default parameter for the other PDFs (e.g. exponential) and how to take into account the hardness of the problem to define this value (e.g. in terms of the landscape ruggedness of the instance to solve). Moreover, we plan to explore portfolio approaches, i.e. ones which combine multiple solvers as well as experiment with techniques for parameter auto-tuning. This experimentation entails a deep analysis of the parallel performance behavior.

2 Fig. 1 :

 21 Fig. 1: CPLS framework structure

Fig. 2 :

 2 Fig. 2: Different probability distribution functions for EO

Table 3 :

 3 Cooperative parallelism on 128 cores (timeout = 300 s)

Source code and instances are available from cri-hpc1.univ-paris1.fr/qap/.

With a probability pAdopt.

Acknowledgments

The authors wish to acknowledge Stefan Boettcher (Emory University) for his explanations about the Extremal Optimization method. The experimentation used the cluster of the University of Évora, which was partly funded by grants ALENT-07-0262-FEDER-001872 and ALENT-07-0262-FEDER-001876.

half of the 66 problem instances now become fully solved (the average time being 3.2s). However, among the 33 remaining ones, 19 remain never solved and, on average, the remaining 14 get solved 6 out of 10 times. Moreover, the average APD over the not fully solved 33 is 0.372%. This contrasts with the previous situation (2.2%), which indicates a significant improvement in the quality of solutions: even when the optimum is not reached, the solution which was found is close.

Stage 3: Cooperative Parallelism

In this final experiment, we attacked the 33 hardest instances with parallelism and cooperation. This was simplified thanks to the CPLS framework which provides the necessary abstraction layers and already handles the communication (see Section 2.3). The sequential EO-QAP needed a very simple adaptation: every R iterations it has to send its current configuration to the Elite Pool and, every U iterations, it retrieves a configuration from the pool, which it nondeterministically 4 adopts it if it is better than the current one. The CPLS system


Danny Munera
email: danny.munera@malix.univ-paris1.fr

Daniel Diaz
email: daniel.diaz@univ-paris1.fr

Salvador Abreu

Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem

Keywords: QAP, heuristics, parallelism, cooperation, hybridization, portfolio

The Quadratic Assignment Problem is at the core of several real-life applications. Finding an optimal assignment is computationally very difficult, for many useful instances. The best results are obtained with hybrid heuristics, which result in complex solvers. We propose an alternate solution where hybridization is obtain by means of parallelism and cooperation between simple single-heuristic solvers. We present experimental evidence that this approach is very efficient and can effectively solve a wide variety of hard problems, often surpassing state-of-the-art systems.

Introduction

The Quadratic Assignment Problem (QAP) was introduced in 1957 by Koopmans and Beckmann [START_REF] Koopmans | Assignment Problems and the Location of Economic Activities[END_REF] as a model of a facilities location problem. This problem consists in assigning a set of n facilities to a set of n specific locations minimizing the cost associated with the flows of items among facilities and the distance between them. This combinatorial optimization problem has many other real-life applications: scheduling, electronic chipset layout and wiring, process communications, turbine runner balancing, data center network topology, to cite but a few [START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF]. This problem is known to be NP-hard and finding effective algorithms to solve it has attracted a lot of research in recent years. To tackle problems of medium or large size (n > 30) one must resort to incomplete methods which are designed to quickly provide good, albeit potentially sub-optimal, solutions. This is the case of metaheuristics. Since the mid-1980s several metaheuristics have been successfully applied to the QAP: tabu search, simulated annealing, genetic algorithms, GRASP, ant-colonies [START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF]. For solving the hardest instances, the current trend is to specialize existing heuristics [START_REF] James | Multistart Tabu Search and Diversification Strategies for the Quadratic Assignment Problem[END_REF][START_REF] Benlic | Breakout Local Search for the Quadratic Assignment Problem[END_REF] often by combining different metaheuristics (hybrid procedures) [START_REF] Drezner | The Extended Concentric Tabu for the Quadratic Assignment Problem[END_REF][START_REF] Drezner | Extensive experiments with hybrid genetic algorithms for the solution of the quadratic assignment problem[END_REF] and/or to resort on parallelism [START_REF] James | A Cooperative Parallel Tabu Search Algorithm for the Quadratic Assignment Problem[END_REF][START_REF] Tosun | On the Performance of Parallel Hybrid Algorithms for the Solution of the Quadratic Assignment Problem[END_REF].

We recently proposed a sequential Extremal Optimization (EO) procedure for QAP which performs well on the QAPLIB instances [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF]. We developed a cooperative parallel version of this method, a process which was eased thanks to our Cooperative Parallel Local Search (CPLS) framework [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF][START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF] for which we developed an implementation in the X10 programming language [START_REF] Charles | X10: An Object-Oriented Approach to Non-Uniform Cluster Computing[END_REF][START_REF] Saraswat | A Brief Introduction to X10 (for the High Performance Programmer)[END_REF]. This solver (called ParEO) behaves very well on the set of 33 hardest and largest instances of QAPLIB. Using 128 cores and within a short time limit of 5 minutes, ParEO is able to find the best known solution (BKS) in each replication for 15 problems. Only for 8 instances is the BKS never reached. Recent research shows that the most promising way to improve QAP resolution is to resort to hybrid procedures, in order to benefit from the strengths of different classes of heuristics. Such is the case of hybrid genetic algorithms (a.k.a memetic algorithms) [START_REF] Drezner | Extensive experiments with hybrid genetic algorithms for the solution of the quadratic assignment problem[END_REF]. The price to pay for this improvement is a significant increase in the complexity of the resulting solver code. In any case, many of the best known existing methods can be easily parallelized thanks to our CPLS framework.

In this paper we propose an alternative approach for hybridization: we resort to cooperation and parallelism to get "the best of both worlds". To this end, the parallel instances of different heuristics communicate their best solutions during execution, and are able to forgo the current computation and adopt a better solution (hoping it will converge faster). The expected behavior is that a solution which appears to be stagnating inside one solver can be improved by another heuristic. When the second solver can no longer improve on this (imported) solution, maybe the original one can, once again, improve the solution yet a bit more, and so on. It is worth noticing that when the first solver sends its current solution, it continues to work on it until it adopts an external solution, itself. This cooperative portfolio approach behaves like a hybrid solver while retaining the original simplicity of each solver. This is particularly true inside the CPLS framework since solvers need not be aware about the (nature of) other solvers.

We implemented such a hybrid solver on top of the X10 version of CPLS, combining two different solvers: Taillard's robust tabu search (RoTS) [START_REF] Taillard | Robust Taboo Search for the Quadratic Assignment Problem[END_REF] and our EO-QAP [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF] method. We have chosen these two solvers because they are simple and also because it turned out that they present complementary strengths: roughly speaking RoTS is stronger in intensifying the search in a given region while EO-QAP is better at widely diversifying the search. The resulting hybrid cooperative solver (called ParEOTS) displays very good performance, as we shall see further. We show that it scales very well, exhibiting a linear speedup when increasing the number of cores. This solver behaves much better than the cooperative versions of both EO-QAP and RoTS alone. On the 33 hardest instances of QAPLIB, using 128 cores and a time limit of 5 minutes, ParEOTS is able to find the BKS for 26 problems at each replication. Even for the 7 other problems, the quality of returned solutions (measured as a percentage of average solution over the BKS), is significantly improved. We also test ParEOTS on Palubeckis' InstXX instances and on Drezners dreXX instances. Moreover we provide optimal solutions for several InstXX instances and for dre90, dre100 and dre132.

The rest of the paper is organized as follows: section 2 discusses QAP, RoTS and EO-QAP. Section 3 presents our parallel hybrid solver. Several experimental results are laid out and discussed in section 4 and we conclude in Section 5.

In this section we recall some background topics: the Quadratic Assignment Problem (QAP) and the two heuristics we plan to combine: RoTS and EO-QAP.

QAP

Since its introduction in 1957, QAP has been widely studied and several surveys are available [START_REF] Burkard | Quadratic Assignment Problems[END_REF][START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Loiola | A survey for the quadratic assignment problem[END_REF][START_REF] Bhati | Quadratic Assignment Problem and its Relevance to the Real World: A Survey[END_REF]. A QAP problem of size n consists of two n × n matrices (a ij ) and (b ij ). Solving such a problem consists in finding a permutation π of {1, 2, . . . n}, minimizing the objective function:

F (π) = n i=1 n j=1 a ij • b πiπj .
In facility location problems, the a matrix represents inter-facility flows and b encodes the inter-location distances. Moreover, QAP can be also used to model scheduling, chip placement and wiring on a circuit board, to design typewriter keyboards, for process communications, for turbine runner balancing among many other applications [START_REF] Commander | A survey of the quadratic assignment problem, with applications[END_REF][START_REF] Zaied | A Survey of Quadratic Assignment Problems[END_REF].

QAP is computationally very difficult: it is a discrete problem, the objective function contains products of variables and the theoretical search space of an instance of size n has a size n!. QAP has been proved to be NP-hard [START_REF] Sahni | P-Complete Approximation Problems[END_REF] (the traveling salesman problem can be formulated as a QAP) and there is noapproximation algorithm for QAP (unless P=NP). In practice, this means that QAP is one of the toughest combinatorial optimization problems, and one with several real-life applications.

QAP can be (optimally) solved with exact methods like dynamic programming, cutting plane techniques and branch & bound algorithms for medium sizes, e.g. n ≤ 30. For larger problems, (meta)heuristics are the most efficient tool. Over the last decades several metaheuristics were successfully applied to QAP: tabu search, simulated annealing, genetic algorithms, GRASP, ant-colonies [START_REF] Said | A Comparative Study of Meta-heuristic Algorithms for Solving Quadratic Assignment Problem[END_REF].

2.2 RoTS : a Tabu Search Procedure for QAP Tabu search as proposed by Glover [START_REF] Glover | Tabu Search[END_REF] has been widely used since the 1990 to tackle QAP. Unquestionably, one of the most important algorithms for QAP is Taillard's robust tabu search [START_REF] Taillard | Robust Taboo Search for the Quadratic Assignment Problem[END_REF] (RoTS). This algorithm uses an adaptive short-term memory for the tabu list by recording the value assigned to a variable for a while (in order to prevent "reverse assignments"). It also uses a clever aspiration criterion (needed to authorize a tabu move to be performed in special circumstances, e.g. if it improves on the best solution found so far). RoTS also incorporates a long-term memory to ensure a form of diversification, by encouraging moves towards to not yet visited regions. RoTS only requires two user parameters to be tuned: the tabu tenure factor (controlling the time an element remains tabu) and the aspiration factor both of which influence the adaptive memory ( [START_REF] Taillard | Comparison of iterative searches for the quadratic assignment problem[END_REF] provides good references values for these parameters). In practice RoTS is tremendously effective on a wide variety of QAP instances, being able to quickly find high quality solutions. Several BKS for QAPLIB instances have been discovered and/or improved by RoTS. A key feature explaining its speed is that the cost of a solution resulting from a swap can be computed incrementally and further optimized using a tabling mechanism. This results in an evaluation in O(n 2 ), while the naïve algorithm is in O(n 3 ). In addition, Taillard put the source code in the public domain. All these reasons explain the fact that RoTS is directly or indirectly at the root of many other methods to solve QAP [START_REF] James | Multistart Tabu Search and Diversification Strategies for the Quadratic Assignment Problem[END_REF][START_REF] Misevicius | A Tabu Search Algorithm for the Quadratic Assignment Problem[END_REF].

2.3 EO-QAP : an Extremal Optimization procedure for QAP Extremal Optimization (EO) is a metaheuristics inspired by self-organizing processes often found in nature [START_REF] Boettcher | Nature's way of optimizing[END_REF][START_REF] Boettcher | Extremal Optimization: an Evolutionary Local-Search Algorithm[END_REF][START_REF] Boettcher | Extremal Optimization[END_REF]. EO is based on the concept of Self-Organized Criticality (SOC) initially proposed by Bak and on the Bak-Sneppen's model [START_REF] Bak | Punctuated equilibrium and criticality in a simple model of evolution[END_REF]. In this model of biological evolution, species have a fitness ∈ [0, 1] (0 representing the worst degree of adaptation). At each iteration, the species with the worst fitness is eliminated (or forced to mutate). This change affects its fitness but also the fitness of all other species connected to this "culprit" element. This results in an extremal process which progressively eliminates the least fit species (or forces them to mutate). Repeating this process eventually leads to a state where all species have a good fitness value, i.e. a SOC. EO follows this line: it inspects the current configuration (assignment of variables), selects one of the worst variables (according to their fitness) to mutate. For this, it ranks the variables in increasing order of fitness (the worst variable has thus a rank k = 1) and then resorts to a Probability Distribution Function (PDF) over the ranks k to chose the culprit element. This PDF introduces uncertainty in the search process. The original EO proposes a power-law:

P (k) = k -τ (1 ≤ k ≤ n).
This PDF takes a single parameter τ which is problem-dependent. Depending on the value of τ , EO provides different search strategies from pure random walk (τ = 0) to deterministic (greedy) search (τ → ∞). With an adequate value for τ , EO cannot be trapped in local minima since any variable is susceptible to mutate (even if the worst are privileged). This parameter can be tuned by the user (a default value is τ = 1 + 1 ln(n) ). EO displays several a priori advantages: it is a simple metaheuristic, it is controlled by only one free parameter (a fine tuning of several parameters becomes quickly tedious) and it does not need to be aware about local minima. Nevertheless, EO has been successfully applied to large-scale optimization problems like graph bi-partitioning, graph coloring or the traveling salesman problem [START_REF] Boettcher | Extremal Optimization: an Evolutionary Local-Search Algorithm[END_REF].

Recently, we proposed EO-QAP: an EO procedure for QAP [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF]. One notable extension we brought to the original EO is to propose different PDFs and to allow the user to chose the most adequate one for a given problem. The sequential procedure performs well on the whole set of QAPLIB instances: 68 instances are solved (i.e. the BKS could be reached) at each execution, 41 are solved at least once and 25 never. The independent parallel version improves the situation significantly: 33 additional instances are systematically solved, 14 are partially solved and 19 remain unsolved. To tackle this remaining set of 33 instances (14+19) we experimented with cooperative parallelism (this version is called ParEO). In the same time limit, ParEO is able to systematically solve 15 new instances and 18 are solved at least once (8 remain unsolved).

We propose an alternative approach for constructing hybrid search methods, resorting on our Cooperative Parallel Local Search Framework (CPLS) [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF][START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF], to provide the hybridization. In a nutshell, the procedure amounts to having several workers, each following its own strategy, some of which are significantly different from others. The cooperative framework oversees every worker, and makes it possible for it to contribute and benefit from the global effort, by managing a pool of best solution candidates (the elite pool ). The fact that the framework is parallel entitles it to obtain performance benefits by just increasing the count of compute units (cores.) Moreover, the workers themselves need to have little or no knowledge of the environment they are running under.

To test these ideas, we experimented with a solver for QAP -an admittedly difficult problem -for which the individual metaheuristic we chose are our EO-QAP algorithm and the RoTS method.

Cooperative Parallel Local Search

Parallel local search methods have been proposed in the past [START_REF] Alba | Parallel Metaheuristics: A New Class of Algorithms[END_REF][START_REF] Alba | Parallel Metaheuristics: Recent Advances and New Trends[END_REF][START_REF] Diaz | Parallel Constraint-Based Local Search on the Cell/BE Multicore Architecture[END_REF]]. Here we focus on multi-walk methods (also called multi-start) which consist in a concurrent exploration of the search space, either independently or cooperatively, the latter being achieved with communication between processes. The Independent Multi-Walks method (IW) [START_REF] Verhoeven | Parallel Local Search[END_REF] is easiest to implement since the solver instances need not communicate with each other. However, the resulting gain tends to flatten when scaling beyond about a hundred processors [START_REF] Caniou | Large-scale parallelism for constraint-based local search: the costas array case study[END_REF], largely because the inherent diversity which brings about the speedups is not sufficient. In the Cooperative Multi-Walks (CW) method [START_REF] Toulouse | Systemic Behavior of Cooperative Search Algorithms[END_REF], the solver instances exchange information (through communication), hoping to hasten the search process. However, the design and implementation of an efficient such method is a very challenging task: choices abound concerning the communication which impact each other, many of which are problem-dependent [START_REF] Toulouse | Systemic Behavior of Cooperative Search Algorithms[END_REF].

We designed the Cooperative Parallel Local Search (CPLS) framework [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF][START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. This framework, available as an open source library in the X10 programming language, allows the programmer to tune the search process through an extensive set of parameters which, at present, statically condition the execution. CPLS augments the IW strategy with a tunable communication mechanism, which allows for the cooperation between the multiple solver instances to seek either an intensification or diversification strategy in the search. At present, the tuning process is done manually: we have not yet experimented with parameter self-adaptation in the CPLS framework (still an experimental feature).

The basic component of CPLS is the explorer node which consists in a local search-based solver instance. The point is to use all the available processing units by mapping each explorer node to a physical core. Explorer nodes are grouped into teams, of size N P T (see Figure 1). This parameter is directly related to the trade-off between intensification and diversification. N P T can take values from 1 to the maximum number of cores. When N P T is equal to 1, the framework coincides with the IW strategy, it is expected that each 1-node team be working on a different region of the search space, without any effort to seek parallel intensification. When N P T is equal to the maximum number of nodes (creating only 1 team in the execution), the framework is mainly geared towards parallel intensification (however a certain amount of diversification is inherently provided by parallelism, between 2 cooperation actions). Each team seeks to intensify the search in the most promising neighborhood found by any of its members. The parameters which guide the intensification are the Report Interval (R) and Update Interval (U ): every R iterations, each explorer node sends its current configuration and the associated cost to its head node (report event). The head node is the team member which collects and processes this information, retaining the best configurations in an Elite Pool (EP ) whose size |EP | is parametric. Every U iterations, explorer nodes randomly retrieve a configuration from the EP , in the head node (update event). An explorer node may adopt the configuration from the EP , if it is "better" than its own current configuration, with a probability pAdopt. Simultaneously, the teams implement a mechanism to cooperatively diversify the search, i.e. they try to extend the search to different regions of the search space.

Typically, each problem benefits from intensification and diversification to some extent. Therefore, the tuning process of the CPLS parameters seeks to provide an appropriate balance between the use of the intensification and diversification mechanisms, in hope of reaching better performance than the noncooperative parallel solvers (i.e. independent multi-walks). A detailed description of this framework may be found in [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF].

Using the CPLS Framework for Hybridization

The current X10 implementation of the CPLS framework already supports the use of multiple metaheuristics. Adding a new one is simple because CPLS provides useful abstraction layers and handles communication. Adding a new metaheuristic comes down to slightly adapt the sequential algorithm: every R iterations it has to send its current configuration to the Elite Pool and, every U iterations, it needs to retrieve a configuration from the pool, which it may subsequently adopt (with probability pAdopt), should it be better than the current one. The overall resulting solver is thus composed of several instances of the same metaheuristic running in parallel, which cooperate by communicating in order to faster converge to a solution. To date, CPLS includes cooperative parallel versions of three different methods: Adaptive Search, Extremal Optimization and Tabu Search. In the present work, we go one step beyond and propose a new usage of the CPLS framework in order to obtain an hybrid parallel solver. For this, individual workers run instances of different metaheuristics, while still collaborating by communicating with the head node. The basic idea of running different metaheuristics in parallel exchanging elite solutions has been proposed [START_REF] Alba | Parallel Metaheuristics: A New Class of Algorithms[END_REF][START_REF] Talukdar | Asynchronous Teams: Cooperation Schemes for Autonomous Agents[END_REF] but only from a general and theoretical point of view. This can also be viewed as a portfolio approach [START_REF] Gomes | Algorithm portfolios[END_REF] augmented with cooperation.

We chose to experiment with this form of hybridization on QAP combining two metaheuristics: our EO-QAP procedure and the RoTS method, resulting in a solver we call ParEOTS. The communication strategies of CPLS remain unchanged, ensuring cooperation between the explorers which now happen to be running different methods. Figure 2 presents possible interactions due to cooperation and the implementation of the hybrid strategy. The team's EP will now contain configurations stemming from explorers running different heuristics. Here is a possible scenario: inside the same team, an instance E 1 of EO-QAP reports a good configuration C 1 to the EP. Later, an instance R 1 of RoTS retrieves C 1 , improves on it (RoTS being strong at intensification) and obtains a better configuration C 2 , on which it reports back to the EP. Later, C 2 gets adopted by an instance E 2 which, being in a diversification phase, moves to a faraway search region, which may provide yet better solutions. Obviously, other scenarios are possible, e.g. when another EO-QAP explorer E 3 also retrieves C 1 (provided by EO) it gives a "second chance" to this configuration (due to its internal stochastic state it can further improve this configuration). The whole system behaves as a hybrid solver, benefiting from cross-fertilization due to the inherent diversity of the search strategies.

Report Event Head

E-1 EO E-2 EO E-3 RTS E-4 RTS RTS3 Elite Pool Head E-1 EO E-2 EO E-3 RTS E-4 RTS

Elite Pool

Update Event

In this section we present an experimental evaluation of our hybrid parallel method (source code, instances and new solutions will be soon available from http://cri-hpc1.univ-paris1.fr/qap/). All experiments have been carried out on a cluster of 16 machines, each with 4 × 16-core AMD Opteron 6376 CPUs running at 2.3 GHz and 128 GB of RAM. The nodes are interconnected with InfiniBand FDR 4× (i.e. 56 GBPS). We had access to 4 nodes and used up to 32 cores per node, i.e. 128 cores. Each problem is executed 10 times stopping as soon as the BKS (which is sometimes the optimum) is found. This execution is done with a short time limit of 5 minutes (in case the BKS is not reached). Such experiments give an interesting information about the quality of solutions quickly obtainable. All times are given either in seconds for small values (as a decimal number) or in a human readable form as mm:ss or hh:mm:ss). The relevant CPLS parameters controlling the cooperation are (as per [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]):

-Team Size (N P T ): we fixed it to N P T = 16. There are thus 8 teams composed of 16 explorer nodes ; 8 running a EO-QAP solver and 8 running RoTS solver. This is constant over all problems. We did not yet experiment with other splits. -Report and Update Interval (R and U ): we manually tuned U and usually fix R = U/2. -Elite Pool (EP ): its size is fixed to 4 for all problems.

-pAdopt: is set to 1. Any solver instance receiving a better configuration than its current one always switches to the new one.

Scalability Analysis

We start this experimental evaluation by analyzing the scalability of ParEOTS. Such an analysis is not easy, because if the BKS cannot be reached, the runtime is only bounded by the timeout used. It is thus necessary to only consider problems that can be systematically solved by the EO sequential solver (to have a reference time using 1 core). We selected two instances of QAPLIB which require the longest sequential time: tai35a solved on average in 42.399s and lipa70a solved in 57.737s. We then ran these problems with ParEOTS, varying the number of cores from 2 to 128. Figure 3 presents the speedup data and curves obtained with our algorithm (using a log-log scale). The Ideal curve corresponds to linear speedup: time is halved when the number of cores is doubled. For both problems the speedup is linear. Using 128 cores, the best speedup is 126, obtained for tai35a whose execution time now only requires 0.336s.

Evaluation on QAPLIB

We here evaluate the performance of our hybrid solver ParEOTS on a set of 33 hard instances of QAPLIB. We selected this set because it is the most difficult set for the independent parallel version of our EO procedure [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF]. In addition to Fig. 3: Speedup profile using the Hybrid CPLS on two QAPLIB instances raw performance, and for validation purpose, we also want to assess the gain obtained with the hybrid version compared the cooperative parallel versions of its two components: ParEO and ParRoTS (also written in X10 within CPLS). For this, all 3 systems are run under the same conditions (see Section 4). Obviously, ParEO runs 128 instances of our EO procedure, ParRoTS runs 128 instances of RoTS while ParEOTS executes 64 instances of EO-QAP and 64 of RoTS.

To measure the hybrid performance we focus on the number of BKS found by each parallel solver. When running 50% of EO and 50% of RoTS we define a low threshold (low) as the average of #BKS found by both parallel solvers. This corresponds to what can be normally expected. Below this value, the hybrid solver is ineffective. Above, it already performs well. Moreover, we define a high threshold (high) as the maximum of the #BKS of both solvers. Above this value, the hybrid solver performs at least as well as the best single solver (a hybrid solver without gain would need twice the number of cores to obtain such a performance). Obviously low and high can be generalized to an hybridization involving more than 2 solvers. For a given problem, executed n times, the performance (hperf ) of the hybrid solver reaching #bks times the BKS is defined as follows:

hperf =                      #bks -low low , if #bks < low #bks -low high -low , if low < #bks < high and low = high 1 + #bks -high n -high , if high ≤ #bks and n = high 1, if high = #bks = n (1)
The performance ranges in [- 1 presents the results. The parameters used for EO are the same as in [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF]. For RoTS we generally use a tabu tenure = 8n and an aspiration = 4n 2 . The table reports, for each solver, the number of times out of 10 runs the BKS was reached (#bks), the Average Percentage Deviation (APD) which is relative deviation percentage computed as follows: 100 × Avg-BKS BKS (where Avg is the average of the 10 found costs), the average execution time (average of the 10 wall times for one instance) and the numbers of adoptions done by the winning explorer (#ad.). The performance value is also reported. The last row presents the averages of each column (or sums for #bks columns).

It is worth noticing that the overall performance of the cooperative parallel version of the 2 base solvers using a short time limit is rather good. Even so, the hybrid solver clearly outperforms them. Focusing on #BKS, it provides high performance (hperf ≥ 1) for 32 instances (only for tai50a does it behave worse than its two components). Moreover, in 4 cases it obtains a hperf = 2 corresponding to cases where it performs much better than both individual solvers (to such an extent that it obtains the perfect score #BKS = 10). It found the BKS at each replication for 26 problems; this is much better than ParEO [START_REF] Taillard | Robust Taboo Search for the Quadratic Assignment Problem[END_REF] and ParRoTS [START_REF] Zaied | A Survey of Quadratic Assignment Problems[END_REF]. In only 4 cases, could ParEOTS not reach the BKS: this number is 8 for ParEO and 6 for ParRoTS. It is worth noticing that even in these 4 cases, the hybridization is still effective since the APD is lower than for its components. For instance, on the very difficult problem tai256c, the hybridization cannot solve the problem but the APD is 0.178 while it is around 0.263 for both components. Another remarkable case is tho150, for which the hybridization is very effective. The average APD is now 0.007 (0.144 for ParEO and 0.019 for ParRoTS). In fact, it turns out that this problem could even be solved once.

The "summary" row reports interesting numbers. All in one, the average APD of ParEOTS is 0.034 which is much better than 0.138 for ParEO and 0.059 for ParRoTS. Regarding execution times, it is a good surprise to see that the increase of quality does not hamper the speed. In fact, with an average execution time of 85s the hybrid solver is faster than ParEO (148s) and ParRoTS (113s).

Testing on Palubeckis' instances

In 2000, Palubeckis proposed a new hard problem generator with known optimum [START_REF] Palubeckis | An Algorithm for Construction of Test Cases for the Quadratic Assignment Problem[END_REF] and provided a set of 10 hard instances called InstXX . Few results have been published about experiments with them. Palubeckis reports the best solutions found by a repeated local search procedure (called multi-start descent or MSD). In [START_REF] Wu | An Ant Colony Optimization Algorithm for Quadratic Assignment Problem[END_REF] the authors propose an Ant Colony Optimization algorithm (QAP-ACO) and test it on these instances (in this work these instances are called paluXX ).

We experimented in the same setting as previously: with 128 cores and a time limit of 5 minutes. Table 2 displays the results for 3 solvers. In addition to the APD we also provide the best cost value found among the 10 runs. Data is taken from the aforementioned articles. We also provide execution times for QAP-ACO for indicative purposes. We set out to construct a hybridized solver by resorting to a parallel and cooperative multi-walk scheme, which relies on the CPLS framework to provide both the cooperation and the parallel or distributed execution. As a testbed for the idea, we chose to tackle the Quadratic Assignment Problem, because it is recognized as a very difficult problem of significant practical interest and also because benchmark instances abound in the literature. For this we designed ParEOTS: a hybrid cooperative parallel solver combining two methods: our Extremal Optimization algorithm and Taillard's robust tabu search. This hybrid solver is much more efficient than any of its two individual base solvers. Regarding QAPLIB, our hybrid solver is able to reach the best known solution (BKS) for all instances except 4. In most cases it is even able to systematically find the BKS at each replication. Even then, for the 4 not fully solved hardest instances (tai80a, tai100a, tai150b and tai256c), the solutions obtained are very close to the BKS. We also tested the solver on other hard instances. The results on Palubeckis' instances are very good: for the first time, ParEOTS optimally solved all instances up to a size n = 100 (prior to this work only optimal solutions for n = 20 were known). We discovered optimal solutions for sizes n = 30..100 and 2 new best obtained solutions for n = 150 and n = 200. Regarding Drezner's instances, the results are even better: we discovered optimal solutions for all instances (including dre90, dre110 and dre132). This is the first time that optimal solutions for these 3 instances are published.

From our experiments, it became clear that: (1) the coding effort for building a hybrid solver is much lower with our approach than for existing hybrid algorithms, and (2) the performance gain over competing approaches is very significant. The latter aspect can be construed as a sort of evolutionary algorithm, one which blends phenotypes rather than genotypes, all under the supervision of the cooperative framework. As to the former, the changes needed to fit the CPLS scheme are minimal and very simple.

We plan to further explore portfolio approaches, combining more than two types of solver as well as experimenting with techniques for parameter autotuning. Another line entails the induction of solver multiplicity by presenting several instances of the same solver, but set up with different parameters.

2 Fig. 1 :

 21 Fig. 1: CPLS framework structure

Fig. 2 :

 2 Fig. 2: Hybridization in CPLS : combining EO-QAP and RoTS

Table 1 :

 1 ParEOTS on QAPLIB and comparison with ParEO and ParRoTS Table

		1, 2]. if hperf < 0 the hybrid solver is ineffective
	on that problem. For hperf ∈ [0, 1) the performance is acceptable and when
	hperf ∈ [1, 2] the performance is very good.
	low	high	very high
	-1	0	1	2

Acknowledgments

The authors wish to thank Prof. E. Taillard for providing the RoTS source code and explanations. The experimentation used the cluster of the University of Évora, which was partly funded by grants ALENT-07-0262-FEDER-001872 and ALENT-07-0262-FEDER-001876.

els19 10 0.000 0.0 2.6 1.00 10 0.000 0.0 0.2 10 0.000 0.0 0.3 kra30a 10 0.000 0.0 3.9 1.00 10 0.000 0.0 2.6 10 0.000 0.0 3.6 sko56 10 0.000 1.5 0.3 1.00 10 0.000 4.8 2.5 10 0.000 0.6 0.0 sko64 10 0.000 1.7 0.3 1.00 10 0.000 4.8 1.5 10 0.000 1.3 0.0 sko72 10 0.000 8.7 1.0 1.00 10 0.000 0:13 1.4 10 0.000 0:16 1.7 sko81 10 0.000 0: [START_REF] Boettcher | Nature's way of optimizing[END_REF] Even with a limit of 5 minutes, the performances of ParEOTS are very good. The optimum is reached for problems whose size n ≤ 100. In addition, for all n ≤ 80 ParEOTS reaches the optimum at each replication. For sizes n > 100, clearly a limit of 5 minutes is too short to reach the optimum. Nevertheless, the obtained solutions are of good quality with an APD around 1.12%: 2-3 times better than challengers. It is worth noticing that for n > 20 all published best obtained solutions are improved (in bold font in the table). Regarding execution times, ParEOTS also outperforms its competitors.

Testing on Drezner's instances

In 2005, Drezner and al. designed new QAP instances with known optimum but specifically ill conditioned to be difficult for metaheuristic methods [START_REF] Drezner | Recent Advances for the Quadratic Assignment Problem with Special Emphasis on Instances that are Difficult for Meta-Heuristic Methods[END_REF]. The authors reports the best solutions found by a powerful compounded hybrid genetic algorithm (called CHG in what follows). The instances are really difficult and only very recently were some results published by Acan and Ünveren with a great deluge algorithm (called TMSGD) [START_REF] Acan | A Great Deluge and Tabu Search Hybrid with Two-stage Memory Support for Quadratic Assignment Problem[END_REF]. These hard instances are thus an interesting challenge for our hybrid solver.

We ran it under the same conditions as before: using 128 cores and with a time limit of 5 seconds. Table 3 presents the results for 3 solvers. Data is taken from the above mentioned articles (in the case of CHG each problem was executed 20 times, presented #BKS are divided by 2 for normalization). We also provide execution times for TMSGD for indicative purposes (TMSGD was run on a 2.1 GHz PC). The performance of ParEOTS is very good: all problems could be optimally solved, and, to the best of our knowledge, this is the first time that an optimal solution is found for dre90, dre110 and dre132. TMSGD performs better than CHG (but the CHG experiment is old). Regarding execution times, ParEOTS needs 2:47 to solve dre90 while TMSGD cannot solve it even using 1:36:33 (CHG reports one hour for dre90 and also fails to find the optimum).


M Slazhniev

A Kizilova

K Bogdan

Kyung Hyun Kim

Hyun Suk Sim

New concept of a modern magnetodynamic installation based on the principle of superposition of pulsating electromagnetic forces

Keywords: electromagnetic forces, superposition, magnetodynamic installation, pressure, modeling

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Known on today existing magnetodynamic installations (MDI), as an induction channel furnace with electromagnetic pump, have a widely application in industrial casting production at melting and complex processing of non-ferrous and ferrous metals and alloys [1,[START_REF] Lee | Proceedings of the 12th International Conference on Aluminium Alloys[END_REF]. Such equipment has been developed in the Physico-Technological Institute of Metals and alloys of National Academy of Sciences of Ukraine (PTIMA NASU). The classic solution of MDI, for example by magnetodynamic installation for aluminum alloys MDN-6A, is designed with dual combined W-shaped channel, with two independently controlled inductors for each channel, and additional electromagnet, for creation of external alternating magnetic field (B). The poles of C-shape electromagnet are located in place connecting of bottom horizontal induction channels with central branches (pumping pipe), forming T-shape working area (WA). The main working principle of this MHD-equipment is based on the interaction of AC electric current (I) with external AC magnetic field (B), to provide to generating of unidirectional pulsating electromagnetic force and pumping of the molten metal from WA, on metal duct, to casting molds as a EM pump [1]. The basic requirement for existent magnetodynamic installations (MDI) and EM pumps there are still in increasing pressure characteristics (p em -electromagnetic pressure), values of mass transferring flow rate (Q) at controlled electromagnetic casting (pouring). Additionally, at applying MDI for processing Al-alloys to need improve the hydrodynamic efficiency of physical methods influencing on melt for implement of the electrophysical melt processing due to action of high densities of AC current, the maximum value of the melt moving velocity, etc. Issue of raising the operational characteristics of MDI can be solved by increasing the hydrothermal and hydraulic efficiency of working area (on example of traditional three-throw WA), where are to create of the electromagnetic pressure. Known, that in existing MDI for aluminium alloy the real coefficient transformation of the electromagnetic force to hydraulic pressure makes no more than 12%. Improving efficiency of existing MHD-pump possible due to decreasing of hydraulic resistance in WA by means optimizing of it configuration and geometry, eliminate the negative impact of electrodynamic effects, optimizing the spatial (in cross section) redistribution of current density j and magnetic field (B) [START_REF] Dubodelov | 9 th International PAMIR Conference on Fundamental and Applied MHD[END_REF]. However, by the theoretical calculation and estimations, such approaches (named below) can provide to increase the magnitudes of electromagnetic force (in WA) on 15 ÷ 25% and a hydraulic pressure up to 30%. But, increasing the amplitude values of the AC current density (j) in the WA is limited by critical values (for molten Al alloys -j < j critical(Al) = 25 A/mm 2 ), that at the exceeding to make the risky appearing of the theta-pinch effect in melt into channel [START_REF] Slazhniev | [END_REF]. Increasing the intensity of the magnetic field in WA requires the development of expensive and large size complicated electromagnets, and can be provide to appearing negative "edge effect" [1]. A radical solution for increasing the efficiency of induction MHD-pump for molten metals, which use AC electrical power for induce of electric current and an external magnetic field, to consist in creation the conditions of vector superposition (addition) of several pulsating unidirectional electromagnetic forces in one local WA of MDI. Due to perform increase the total component of electromagnetic force, which create the hydraulic (static) pressure in molten melt (static), at simultaneously decreasing useless action of vibration (pulsation) component.

Presentation of the problem

According to the canons of the MDI operating theory [1] the electromagnetic pressure there is as action result of the created in MDI of volumetric electromagnetic forces (EMF) (1) to output branch of three-throw working area WA (or cruciform working area for some known solutions MHD-pumps) [1]: [START_REF] Lee | Proceedings of the 12th International Conference on Aluminium Alloys[END_REF] where j m и B mpeak values of j (density of electric AC), A/mm 2 , and B (magnetic field), Т; ω -angular frequency rad/sec; φangle of shift phase between j and B, rad; p em -electromagnetic pressure (2), Pa; a, b -the length and width WA, m; Re -real part of the vector product; Vvolume of WA, m 3 ; l -height of WA, m.

                     t 2 cos B j 5 . 0 cos B j 5 . 0 t sin B ) t sin( j B j dV df F m m m m m m em em.f. (1)       cos jBl 2 dV * jB Re ab 4 1 p V em
In equations ( 1) the constant component 0.5j m B m cos(φ) to create the pressure but component 0.5j m B m cos(2ωt + φ), provides of doubled frequency vibratory action in melt (2ω = 2πf). Thus, in the existing MDI for making hydraulic pressure in melt has been use only 50% of peak value of generated in WA MDI pulsating EMF [1]. By proposed in this study "concept" the improving efficiency of pulsing electromagnetic force transformation to hydraulic pressure in MDI (as in electromagnetic pumps for molten melts) there is possible at making superposition (vector summation) of three vertical directed pulsating electromagnetic forces (F em.f.1 , F em.f.2 , F em.f.3 ), which have phase shift relative to each other  phase = 2π. At them each pulsing EMF with own angle of phase shift can be create by own localized single or two-phase alternating currents AC (I) and magnetic field (B), which also have a phase shift relative to each other  phase = 2π/3. Thus, the resulting superposition of three vertically unidirectional pulsating electromagnetic forces (PEMF) in one local WA provide to create in three time higher component of EMF, which make a hydraulic pressure in melt (∑F em.f. = 1.5j m B m cos(φ)), but vibrating useless action, due to phase shift angle to each other (2 I = 0, 4π/3, 8π/3), have a minimal values and cancel each (Fig. 1). For cases that I 1 = I 2 = I 3 = I m and B 1 = B 2 = B 3 = B m there is vibration component equal "0". Fig. 1: Graph result of superposition (∑F em.f. ) in three-throw one local working area (WA) of three unidirectional pulsating electromagnetic forces (F em.f.1 , F em.f.2 , F em.f.3 ), with phase shift between each one -120 ( 1 = 0,  2 = 2π/3;  3 = 4π/3).

) t 2 cos( B I 5 , 0 ) cos( B I 5 , 0 ) t cos( B ( )) t cos( I ( F 1 1 1 1 1 1 em.f.1                    ; ) 3 4 t 2 cos( B I 5 , 0 ) cos( B I 5 , 0 ) 3 2 t cos( B ( )) 3 2 t cos( I ( F 2 2 2 2 2 2 em.f.2                          ; (3) ) 3 8 t 2 cos( B I 5 , 0 ) cos( B I 5 , 0 ) 3 4 t cos( B ( )) 3 4 t cos( I ( F 3 3 3 3 3 3 em.f.3                            0 ) cos( B I 3 5 . 0 F F F F m m em.f.3 em.f.2 em.f.1 em.f.            (4) F em.f
New original design for the three phase's of magnetodynamic pump and furnace was developed to the technical implementation and experimental verification of the superposition principle (Fig. 2).

As opposed to existing MDI with W-shaped two induction channels, in the proposed here the new design and solution MDI have an induction channel with three vertically connected induction channels, which have one vertical pumping pipe (Fig. 2 and3). Each from three induction channels to connected in the vertical plane on 120 о angle each other (Fig. 2, 3). Horizontal triple channel (Fig. 3) to consist of three horizontal slopes (Fig. 3), which connected at an angle of 90 о with a central/vertical channel and with 120 о angle relative to each other in a horizontal plane. In this case under central pumping channel/pipe to forming the three-throw one local working area (WA).

In addition, as opposed to existing, in new design MDI the three phase's electromagnet (Fig. 3) has three poles, that there are directed at angle of 90 о to three-throw WA and arranged horizontally from three sides at angle of 60 о -relative to lateral horizontal channels and 120 о angle relative to each other. According to the ideas in formed in three-throw WA (Fig. 3) at the interaction of mutually perpendicular horizontal electric currents I* 1 , I* 2 , I* 3 with three alternating magnetic fields (В 1 , В 2 , В 3 ), in pole gap of three poles electromagnet, at switching the windings in three-phase electricity network by a "triangle" -formed three pulsating electromagnetic forces directed vertically (F em.f.1 , F em.f.2 , F em.f.3 ). The AC  1 ( 1 * =  1 +  2 ) has induced by the first and second inductors,  2 ( 2 * =  2 +  3 by second and third inductors and  3 ( 3 * =  3 +  1 ) by first and third inductors at switching of each coils to three-phase electrical network by a "star". Each one electromagnetic force there are to changing in dynamic mode (Fig. 4) in time intervals -t i = π/2 +  i , equal to 1/3 period (T i = 2π/3 =  i ). Next, the process generation and changing for each electromagnetic force (directed vertically) continues in a time cycle -2π, 4π, etc. The model/simulator as a solid coil was developed for experimental verification of proposed idea (Fig. 5). Such model to consist of three horizontal slopes, which were made from aluminum plate 10 mm thick, height up to 80 mm, connected in the horizontal plane at an angle 120). All three horizontal slopes of physical simulation model was connected at 90 angles to central cylindrical aluminum rod (diameter 40 mm), forming in result three-throw WA MDI (Fig. 5).

3-phase's electromagnet of simulation model was performed from three not connected each other horizontal magnetic cores (cross section 70×90 mm and a length of 400 mm), with the single layer low voltage windings (30 turns), and withal are used as the inductors.

For induce of AC in model used induction method at placing by two turns (coil) of copper cable (cross section 120 mm 2 ) on each inductor. The ends (clamps) each coper cable has connected to the ends of the horizontal slopes of model, forming shorted secondary coils in each of the three electric circuits by "star" (Fig. 5). Powering of the electromagnet/inductors in the model the low-voltage transformers from MDN-6A and voltage from 20 to 80 Volts used.

For experimental measuring and recording the electromagnetic force impact, the central vertical rod of solid model/simulator was attached to strain gauges sensor, which measured the weight equivalent of ponderomotive force action (F = mg). At the experimental test of solid simulation model has been performed at creation in the poles gap of electromagnet magnetic field in range from 0.01 to 0.15 T, and the AC current in the horizontal coil respective to own magnet pole 0.

Conclusion

Theoretical, analytical and experimental modeling (using physical and computer simulation) showed efficiency of the developed method of superposition of three pulsing electromagnetic forces in the original developed design of threethrow working area prototype of new design magnetodynamic installations for molten metals for heating, stirring, holding, processing and electromagnetic pumping, dosage and low pressure casting technologies. Shown, that in comparison with existing MDI for aluminium alloys (MDN-6A), by superposition principle to possible increase the hydraulic pressure in molten metals in 3 times from 30 kPa to 100 kPa, at same electric current density (j) in inductions channel and intensity of magnetic field (B) into electromagnet poles gap. New principle of superposition to provide the maximally using of pulsation electromagnetic forces power for transformation its energy to moving molten metals in AC EM pump due to optimal assuming useful and canceling useless components of pulsating EM forces. Obtained technical and theoretical results can be applied at creation new effectively design of modern magnetodynamic installations for non-ferrous and ferrous metals and alloys in modern foundry and casting industries.

2 Fig. 2 :Fig. 3 :

 223 Fig. 2: Scheme of the three phase's prototype of MDI implementing the principle of superposition of PEMF

3 Fig. 4 :

 34 Fig. 4: Dynamic scheme (in time) interleaving of vertical pulsing EMF in WA (t i = π/2 +  i ) by each 1/3 period.

4 to 2

 2 Fig. 5: Experimental model of the triplex inductions channel with three-throw WA (aluminum plate)


A Timofeev

Maksimov

V Demidovich
email: vbdemidovich@mail.ru

M Khatsayuk

V Timofeev

A Maksimov

I Rastvorova

Numerical Simulation of MHD Processes in the Technology of Non-crucible Induction Melting of Titanium Alloys

Keywords: non-crucible melting, titanium alloy, MHD, numerical modelling

come   L'archive ouverte pluridisciplinaire

Introduction

A unique opportunity to obtain a liquid phase of a titanium alloy inside cylindrical bars with induction heating reveals prospects for development of ultimately new technology for non-vacuum melting of titanium alloys. Such a technology looks competitive and energy-efficient compared to the existing technology of induction melting in a cold crucible due to the fact that no additional equipment is used for formation of vacuum and the process of obtaining titanium melt inside of the ingot by induction method requires much less time and energy consumption. Fig. 1 shows the dynamics of temperature field changes along the bar cross-section in the process of induction heating and formation of the melt inside of the ingot. Before the moment of time t 1 there is intensive heating of the bar surface with the constant supplying power value, the surface temperature significantly exceeds the core temperature (T SURFCE > T CORE ). Further the power value decreases, due to thermal losses from the bar surface the temperature maximum starts to shift to some depth from the surface and the temperature drop between the surface and the core starts to even. At the moment of time t 2 the temperatures on the surface and core are equal (T SURFCE = T CORE ) and consequently ΔТ 2 = 0. This effect takes place with induction heating of all metals however it shows itself more for titanium alloys due to low thermal conductivity and high melting temperature. Overheating of inner layers of metal can finally lead to the beginning of their melting at the moment of time t 3 . The melting process takes place before there appears thermodynamic balance between energy coming into the bar and thermal losses from its surface when there is a melt zone inside of the bar separated from the external environment by a layer of protective skull ΔХ (moment of time t 4 ). Implementation of the offered technology requires careful choice of geometric and energetic parameters of the system "inductor -bar" as well as Fig. 1: Temperature field changes dynamics around the bar radius. the parameters of heating modes and hold. Preliminary research of the system was carried out on a thermal model [START_REF] Demidovich | Technologies to electromagnetically produce cast products of titanium and titanium alloys[END_REF]. The results of theoretical research on the present model showed just a partial coincidence with the experiments (Fig. 2). There was a theory that in order to get a correct picture of the process it is necessary to complete the model taking into account magneto hydrodynamic (MHD) effects. At the initial stage it is enough to obtain a solution on a two-dimensional model.

Mathematic model

The system "inductor -bar" has an axial symmetry and the electromagnetic task can be solved in an axis -symmetrical 2D setting. Such systems usually have low values of the magnetic Reynolds number which helps to solve the electromagnetic task not taking into account the influence of velocity fields in the metal on an electromagnetic field. Thus to solve the electromagnetic task the following allowances have been made: electromagnetic field in a calculation area changes by the harmonic law and has an axial symmetry in a 2D setting, we neglect the influence of metal motion on electromagnetic field. The system of equations describing the electromagnetic field consists of equations for vector and scalar potentials and equations of continuity. A further stage of mathematic modelling is solving a thermal hydrodynamic task based on the results obtained during electromagnetic calculation. The system of equations describing thermal hydrodynamic processes consists of the equation of mass, motion and energy conservation. As we should anticipate appearance of free convection currents and their turbulization, it is necessary to introduce gravitation forces taking into account the difference of densities and application of turbulence models. For the account of free convection motions the Boussinesq-Oberbeck approximation is used. In this approximation the dependence of density on temperature gets linear and is only considered with volume forces. Thus, the liquid can be considered incompressible. To take account of the processes taking place during melting and freezing, the thermal hydrodynamic model is completed with a model of crystallization and melting based on the method "enthalpy -porosity". The influence of the electromagnetic field is taken into account by introduction of the corresponding source members into the equations of motion and energy which have been obtained during the electromagnetic calculation. The thermal hydrodynamic task can be solved in an axis symmetrical 2D setting. Today there are several approaches to solve a multi-phase MHD task and several ways of their implementation [START_REF] Baake | Analysis of transient heat and mass transfer processes in the melt of induction channel furnaces using LES[END_REF][START_REF] Pesteanu | The Multicell Volume of Fluid (MC-VOF) Method for the Free Surface Simulation of MFD Flows. Part II: Experimental Verifications and Results[END_REF][START_REF] Bojarevics | Magnetic levitation of large liquid volume[END_REF][START_REF] Spitans | Numerical modelling of free surface dynamics of conductive melt in the induction crucible furnace[END_REF]. An alternative algorithm of solving such tasks has been suggested on the basis of software products ANSYS and Fluent [START_REF] Yu | Induction installation with MHD influence on high aluminum alloy during preparing and casting processes[END_REF][START_REF] Khatsayuk | Mathematical modeling of casting processes in electromagnetic field[END_REF]. Their peculiarity is transfer of full contribution of the demanded value in the calculation area of the thermal hydrodynamic part of the task. Such an approach can be applied for a broad spectrum of MHD tasks including distribution of various phases, the boundaries of which are located in the concentration area of an electromagnetic field and lead to its responsive distortion. Numerical calculation of the thermal hydrodynamic part of the task taking into account the phase transit and turbulent effects has been fulfilled in Fluent, an electromagnetic part in ANSYS Classic. To provide the transfer of source members of the motion and energy equations obtained as a result of solution of the electromagnetic part of the task from ANSYS into Fluent and transfer of the liquid phase distribution, obtained as a result of solving the thermal hydrodynamic part of the task from Fluent into ANSYS, data exchange algorithm has been developed and fulfilled. To fulfill the algorithm the following languages were used: a functional language Scheme, a structural language C with the library User-defined function (UDF) based on the in-built compiler Fluent and a special Ansys parametric design language (APDL). Copper inductor Bar made of alloy ВТ6 Fig. 3: Geometry and the adopted system dimensions.

Research on the dynamics of liquid phase creation

Based on the created mathematic modelling there has been a numeric experiment carried out on formation of a liquid phase inside of a cylindrical bar of titanium alloy ВТ6 in an electromagnetic field of a 10-wind inductor with the current 2500 А with the frequency of 4 kHz. The emissivity factor on the surface of the bar was taken as 0.5. The main geometrical dimensions of the system are shown in Fig. 3.

To simulate the maintenance of the constant surface temperature in the model there is a control algorithm for supplied power. When the surface is heated above the temperature of 1407 °С there is a smooth decrease of supplied power by 30% to the temperature of 1644 °С. The function which controls the power is as follows:

(

where P 0base power; k fall = 0.3coefficient of power decrease; T serf maxmaximum temperature on the surface; T S , T F temperature of the beginning and end of power decrease. The obtained distribution of the liquid phase concentration in the bar at different moments of time is shown in Fig. 4-a.

where Vthe bar volume.

As it is seen from the pictures, the liquid phase formation takes place at a distance from the surface with the following distribution along the whole depth of the bar. The character of the liquid phase formation process at an initial melting stage is defined by distribution of enthalpy and the following mass transfer process. A condition necessary for the beginning of melt, when the temperature of solidus T SOL is reached, is sufficient power reserve to increase specific enthalpy by the value of metal L latent heat. The present power reserve is defined by the difference between active power, emitted in the bar due to eddy currents running on it, and thermal losses from the bar surfaces which are mainly caused by thermal emissions. (Fig. 5).

As electromagnetic efforts influence the bar, there are hydrodynamic currents appearing in the forming liquid phase. Heat and mass transfer taking place by means of those currents defines further dynamics of the liquid phase distribution. The distribution of the temperature field, velocity vector fields, specific electromagnetic force and power at different moments of time are shown in Fig. 6.

As it is seen from the pictures, the currents character is a symmetrical two-contour toroidal circulation typical of classical systems "inductor -melt". As the liquid phase volume increases, the size and intensity of circulation increases. At once the square of the phase transition wetted border increases which speeds up melting of metal central

Conclusion

A coupled numerical model for non-vacuum melting process of titanium alloy ВТ6 has been made. The model takes into account MHD-effects and processes of phase transition with correction of the electromagnetic part of the task. The comparisons of the results obtained on the present model with the results obtained from a thermal electrical model has shown how important it is to consider latent melting heat and forced convection currents in a liquid phase due to the action of an electromagnetic field. These results showed better agreement with experiments presented in [START_REF] Maslikov | Research of obtaining liquid phase of titanium alloys in cylindrical bodies during induction heating[END_REF]. Further, this model will be used to find out key parameters and patterns for transformation of energy in the system. However, for a more detailed research of non-stationary turbulent MHD-effects and their influence on the process of liquid phase formation it is necessary to build a 3D model with the application of LES modelling of turbulence.

Fig. 2 :

 2 Fig. 2: Results of electro-thermal calculation (left) and experiment (right).

Fig. 4 -

 4 b shows the dynamics of the liquid phase percentage in the bar volume. The liquid phase percentage is defined from the equation: ∫

Fig. 4 :

 4 Dynamics of the liquid phase formation.

Fig. 5 :

 5 Fig. 5: Distribution dynamics of integral effective power and power losses in a bar.

  bar surface layers retain their shape and a solid phase thickness providing its tightness and not allowing the melt to leak outside.

Processes in the Technology of Non-crucible Induction Melting of Titanium Alloys. 8th International Conference on Electromagnetic Processing of Materials (EPM 2015), Oct 2015, Cannes, France. hal-

Acknowledgment

Financial support from Russian foundation for basic research under contract 15-38-50067 "Mathematical modeling and research of physical processes during non-crucible melting of titanium alloys in electromagnetic field" is gratefully acknowledged.