M Mezmaz

N Melab

Y Kessaci
email: yacine.kessaci@lifl.fr

Y C Lee
email: yclee@it.usyd.edu.au

E-G Talbi

A Y Zomaya
email: zomaya@it.usyd.edu.au

D Tuyttens
email: daniel.tuyttens@umons.ac.be

A Parallel Bi-objective Hybrid Metaheuristic for Energy-aware Scheduling for Cloud Computing Systems

Keywords: Energy-aware Scheduling, Cloud Computing, Metaheuristics

In this paper, we investigate the problem of scheduling precedence-constrained parallel applications on heterogeneous computing systems (HCSs) like cloud computing infrastructures. This kind of applications was studied and used in many research works. Most of these works propose algorithms to minimize the completion time (makespan) without paying much attention to energy consumption.

We propose a new parallel bi-objective hybrid genetic algorithm that takes into account, not only makespan, but also energy consumption. We particularly focus on the island parallel model and the multi-start parallel model. Our new method is based on dynamic voltage scaling (DVS) to minimize energy consumption.

In terms of energy consumption, the obtained results show that our approach outperforms previous scheduling methods by a significant margin. In terms of completion time, the obtained schedules are also shorter than those of other algorithms. Furthermore, our study demonstrates the potential of DVS.

Hybridization; Parallelization; Genetic algorithm; Precedence-constrained parallel applications

Introduction

Precedence-constrained parallel applications are one of the most typical application model used in scientific and engineering fields. Such applications can be deployed on homogeneous or heterogeneous systems (HCSs) like cloud computing infrastructures.

Cloud computing is a simple concept that has emerged from heterogeneous distributed computing, grid computing, utility computing, and autonomic computing. In cloud computing, end-users do not own any part of the infrastructure. The end-users simply use the services available through the cloud computing paradigm and pay for the used services. The cloud computing paradigm can offer any conceivable form of services, such as computational resources for high performance computing applications, web services, social networking, and telecommunications services. Several criteria determine the quality of the provided service, the production cost of this service, and therefore the price paid by the end-user. The duration of this service (makespan) and the consumed energy are among of these criteria. The idea is to provide end-users a more flexible service that takes into account the duration of the service and the consumed energy. End-users could then find the right compromise between these two conflicting objectives to solve their precedence-constrained parallel applications.

The problem of finding the right compromise between the resolution time and the energy consumed of a precedence-constrained parallel application is a bi-objective optimization problem. The solution to this problem is a set of Pareto points. Pareto solutions are those for which improvement in one objective can only occur with the worsening of at least one other objective. Thus, instead of a unique solution to the problem, the solution to a biobjective problem is a (possibly infinite) set of Pareto points. To the best of our knowledge, there is no research published in the literature to solve the above problem with a Pareto approach.

However, many works have focused on precedence-constrained parallel applications (e.g., [START_REF] Darbha | Optimal scheduling algorithm for distributed-memory machines[END_REF], [START_REF] Zomaya | Genetic scheduling for parallel processor systems: comparative studies and performance issues[END_REF], [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] and [START_REF] Lee | A novel state transition method for metaheuristic-based scheduling in heterogeneous computing systems[END_REF]). Most of these works propose algorithms to minimize the makespan. Only recently that some works are interested in minimizing the energy consumption (e.g., [START_REF] Kim | Power aware scheduling of bag-oftasks applications with deadline constraints on DVS-enabled clusters[END_REF], [START_REF] Bunde | Power-aware scheduling for makespan and flow[END_REF], [START_REF] Zhu | Scheduling with dynamic voltage/speed adjustment using slack reclamation in multiprocessor realtime systems[END_REF], [START_REF] Zhu | Power-aware scheduling for and/or graphs in real-time systems[END_REF], [START_REF] Ge | Performance-constrained distributed dvs scheduling for scientific applications on power-aware clusters[END_REF] and [START_REF] Rountree | Bounding energy consumption in large-scale MPI programs[END_REF]).

Using green IT techniques can significantly reduce an organization's and ultimately a country's carbon footprint. The UN bought 461 tons of carbon offsets to ensure that the September 2009 Summit on Climate Change in New York was carbon neutral. According to [START_REF] Cameron | Trading in green it[END_REF], if green IT techniques were implemented that reduced the energy use of 10,000 computers in West Virginia by 50%, the deployment would produce 33 times less carbon in one year than that produced by the summit. A recent study on power consumption by servers [START_REF] Koomey | Estimating total power consumption by servers in the US and the world[END_REF] shows that, in 2005, the power used by servers represented about 0.6% of total U.S. electricity consumption. That number grows to 1.2% when cooling and auxiliary infrastructures are included. In the same year, the aggregate electricity bill for operating those servers and associated infrastructure was about $2.7 billions and $7.2 billions for the U.S. and the world, respectively. The total electricity use for servers doubled over the period 2000 to 2005 in worldwide. The number of transistors integrated into today's Intel Itanium 2 processor reaches nearly 1 billion. If this rate continues, the heat (per square centimeter) produced by future Intel processors would exceed that of the surface of the sun [START_REF] Koch | Discovering multi-core: Extending the benefits of Moore's law[END_REF]. This implies the possibility of worsening system reliability, eventually resulting in poor system performance.

Due to the importance of energy consumption, various techniques including dynamic voltage scaling (DVS), resource hibernation, and memory optimizations have been investigated and developed [START_REF] Venkatachalam | Power reduction techniques for microprocessor systems[END_REF]. DVS among these has been proven to be a very promising technique with its demonstrated capability for energy savings (e.g., [START_REF] Bunde | Power-aware scheduling for makespan and flow[END_REF], [START_REF] Zhu | Scheduling with dynamic voltage/speed adjustment using slack reclamation in multiprocessor realtime systems[END_REF] and [START_REF] Ge | Performance-constrained distributed dvs scheduling for scientific applications on power-aware clusters[END_REF]). For this reason, we adopt this technique and it is of particular interest to this study. DVS enables processors to dynamically adjust voltage supply levels (VSLs) aiming to reduce power consumption. However, this reduction is achieved at the expense of sacrificing clock frequencies.

In this paper, we investigate the energy issue in task scheduling particularly on HCSs like cloud computing systems. We propose a new parallel bi-objective hybrid genetic algorithm that takes into account, not only makespan, but also energy consumption. Our new approach is a hybrid between a multi-objective parallel genetic algorithm and energy-conscious scheduling heuristic (ECS) [START_REF] Lee | Minimizing energy consumption for precedence-constrained applications using dynamic voltage scaling[END_REF]. The results clearly demonstrate the superior performance of ECS over the other algorithms like DBUS [START_REF] Bozdag | A task duplication based bottom-up scheduling algorithm for heterogeneous environments[END_REF] and HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF]. Genetic algorithms make it possible to explore a great range of potential solutions to a problem. The exploration capability of the genetic algorithm and the intensification power of ECS are complementary. A skillful com-bination of a metaheuristic with concepts originating from other types of algorithms lead to more efficient behavior.

Our algorithm is effective as it profits from the exploration power of the genetic algorithm, the intensification capability of ECS, the cooperative approach of the island model, and the parallelism of the multi-start model. The island model and the hybridization improve the quality of the obtained results. The multi-start model reduces the running time of a resolution. Furthermore, one of the major interests of our approach is to give the enduser a set of Pareto solutions to choose according to the desired quality of service, in particular the completion time, and the cost of the service in terms of energy and consequently in terms of price willing to pay. The proposed method can easily be applied to loosely coupled HCSs (e.g., cloud computing systems) using advance reservation and various sets of frequency-voltage pairs.

Our new approach is evaluated with the Fast Fourier Transformation task graph which is a real-world application. Experiments show that (1) the hybridization improves on average the best known results obtained in the literature (by 47.5% for the energy consumption and 12% for the completion time), (2) the island model significantly improves the results obtained using only the hybridization, and (3) the multi-start model accelerates our approach with an average speedup of 13 using 21 cores.

The remainder of the paper is organized as follows. Section 2 presents the application, system, energy and scheduling models used in this paper. Section 3 describes the related work. Our algorithm is presented in Section 4. The results of our comparative experimental study are discussed in Section 5. The conclusion is drawn in Section 6. The paper ends with an appendix which describes our approaches using pseudo-code.

Problem modeling

In this section, we describe the system, application, energy and scheduling models used in our study.

Cloud computing model

A cloud computing system is a set of resources designed to be allocated ad hoc to run applications. In our model, the cloud is assumed to be hosted in a data center which is composed by heterogeneous machines. This data center provides a set of services hosted on thousands of high-end computing servers.

The need in terms of services of an application can be modeled by a task graph. In this graph, an edge between two tasks represents an inter-service communication.

The cloud computing system used in this work consists of a set P of p heterogeneous processors/machines. Each processor p j ∈ P is DVS-enabled; in other words, it can operate with different VSLs (i.e., different clock frequencies). For each processor p j ∈ P , a set V j of v VSLs is random and uniformly distributed among three different sets of VSLs (Table 1). Since clock frequency transition overheads take a negligible amount of time (e.g., 10µs-150µs [11], [START_REF] Min | Dynamic voltage scaling techniques for distributed microsensor networks[END_REF]), these overheads are not considered in our study. The inter-processor communications are assumed to perform with the same speed on all links without contentions. It is also assumed that a message can be transmitted from one processor to another while a task is being executed on the recipient processor which is possible in many systems.

Application model

Parallel programs can be generally represented by a directed acyclic graph (DAG). A DAG, G = (N, E), consists of a set N of n nodes and a set E of e edges. A DAG is also called a task graph or macro-dataflow graph. In general, the nodes represent tasks partitioned from an application; the edges represent precedence constraints. An edge (i, j) ∈ E between task n i and task n j also represents inter-task communication. A task with no predecessors is called an entry task, n entry , whereas an exit task, n exit , is one that does not have any successors. Among the predecessors of a task n i , the predecessor which completes the communication at the latest time is called the most influential parent (MIP) of the task denoted as MIP (n i ). The longest path of a task graph is the critical path.

The weight on a task n i denoted as w i represents the computation cost of the task. In addition, the computation cost of the task on a processor p j , is denoted as w i,j and its average computation cost is denoted as wi .

The weight on an edge, denoted as c i,j represents the communication cost between two tasks, n i and n j . However, a communication cost is only required when two tasks are assigned to different processors. In other words, the communication cost when tasks are assigned to the same processor is zero and thus can be ignored.

The earliest start time of, and the earliest finish time of, a task n i on a processor p j is defined as

EST (n i , p j ) = 0 if n i = nentry EF T (M IP (n i ), p k ) + c M IP (n i ),i otherwise EF T (n i , p j ) = EST (n i , p j ) + w i,j
Note that the actual start and finish times of a task n i on a processor p j , denoted as AST (n i , p j ) and AF T (n i , p j ) can be different from its earliest start and finish times, EST (n i , p j ) and EF T (n i , p j ), if the actual finish time of another task scheduled on the same processor is later than EST (n i , p j ). In the case of adopting task insertion the task can be scheduled in the idle time slot between two consecutive tasks already assigned to the processor as long as no violation of precedence constraints is made. This insertion scheme would contribute in particular to increasing processor utilization for a communication intensive task graph with fine-grain tasks.

A simple task graph is shown in Figure 1 with its details in Table 2 and Table 3. The values presented in Table 3 are computed using two frequently used task prioritization methods, t-level and b-level. Note that, both computation and communication costs are averaged over all nodes and links. The t-level of a task is defined as the summation of the computation and communication costs along the longest path of the node from the entry task in the task graph. The task itself is excluded from the computation. In contrast, the b-level of a task is computed by adding the computation and communication costs along the longest path of the task from the exit task in the task graph (including the task). The b-level is used in this study.

The communication to computation ratio (CCR) is a measure that indicates whether a task graph is communication intensive, computation intensive or moderate. For a given task graph, it is computed by the average communication cost divided by the average computation cost on a target system.

Energy model

Our energy model is derived from the power consumption model in complementary metal-oxide semiconductor (CMOS) logic circuits. The power consumption of a CMOS-based microprocessor is defined to be the summation of capacitive, short-circuit and leakage power. The capacitive power (dynamic power dissipation) is the most significant factor of the power consumption. The capacitive power (P c ) is defined as

P c = ACV 2 f, (1)
where A is the number of switches per clock cycle, C is the total capacitance load, V is the supply voltage, and f is the frequency. Equation ( 1) clearly indicates that the supply voltage is the dominant factor; therefore, its reduction would be most influential to lower power consumption. The energy consumption of the execution of a precedence-constrained parallel application used in this study is defined as

E = n i=0 ACV 2 i f.w * i = n i=0 αV 2 i w * i ,
where V i is the supply voltage of the processor on which task n i is executed, and w * i is the computation cost of task n i (the amount of time taken for n i 's execution) on the scheduled processor.

Scheduling model

The task scheduling problem in this study is the process of allocating a set N of n tasks to a set P of p processors (without violating precedence constraints) that minimizes makespan with energy consumption as low as possible. The makespan is defined as M = max{AF T (n exit )} after the scheduling of n tasks in a task graph G is completed. Although the minimization of makespan is crucial, tasks of a DAG in our study are not associated with deadlines as in real-time systems.

Related work

In this section, we present some noteworthy works in task scheduling, particularly for HCSs, and then scheduling algorithms with power/energy consciousness.

Scheduling in HCSs

Due to the NP-hard nature of the task scheduling problem in general cases [START_REF] Garey | Computers and intractability: A guide to the theory of NP-completeness[END_REF], heuristics, in particular meta-heuristics, are the most popularly adopted scheduling approaches. List scheduling heuristics are the dominant heuristic technique. This is because empirically, list scheduling algorithms tend to produce competitive solutions with lower time complexity compared to algorithms in the other categories [START_REF] Kwok | Benchmarking the task graph scheduling algorithms[END_REF].

The HEFT algorithm [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] is highly competitive in that it generates a schedule length comparable to other scheduling algorithms, with a low time complexity (O(nlogn+(e+n)p)). It is a list-scheduling heuristic consisting of the two typical phases of list scheduling (i.e., task prioritization and processor selection) with task insertion.

The DBUS algorithm [START_REF] Bozdag | A task duplication based bottom-up scheduling algorithm for heterogeneous environments[END_REF] is a duplication-based scheduling heuristic that first performs a critical path based listing for tasks and schedules them with both task duplication and insertion. The experimental results in [START_REF] Bozdag | A task duplication based bottom-up scheduling algorithm for heterogeneous environments[END_REF] show its attractive performance, especially for communication-intensive task graphs. The time complexity of DBUS is in the order of O(n 2 p 2 ).

Scheduling with energy consciousness

Most of previous studies on scheduling with the consideration of energy consumption are conducted on homogeneous computing systems [START_REF] Kim | Power aware scheduling of bag-oftasks applications with deadline constraints on DVS-enabled clusters[END_REF], [START_REF] Zhu | Scheduling with dynamic voltage/speed adjustment using slack reclamation in multiprocessor realtime systems[END_REF], [START_REF] Zhu | Power-aware scheduling for and/or graphs in real-time systems[END_REF], [START_REF] Ge | Performance-constrained distributed dvs scheduling for scientific applications on power-aware clusters[END_REF], [START_REF] Rountree | Bounding energy consumption in large-scale MPI programs[END_REF], [START_REF] Chen | Multiprocessor energy-efficient scheduling for real-time tasks with different power characteristics[END_REF] or single-processor systems [START_REF] Zhong | Energy-aware modeling and scheduling for dynamic voltage scaling with statistical real-time guarantee[END_REF]. In addition to system homogeneity, tasks are generally homogeneous and independent. Slack management/reclamation is a frequently adopted technique with DVS.

In [START_REF] Zhu | Scheduling with dynamic voltage/speed adjustment using slack reclamation in multiprocessor realtime systems[END_REF], several different scheduling algorithms using the concept of slack sharing among DVS-enabled processors were proposed. The rationale behind the algorithms is to utilize idle (slack) time slots of processors lowering supply voltage (frequency/speed). This technique is known as slack reclamation. These slack time slots occur, due to earlier completion (than worst case execution time) and/or dependencies of tasks. The work in [START_REF] Zhu | Scheduling with dynamic voltage/speed adjustment using slack reclamation in multiprocessor realtime systems[END_REF] has been extended in [START_REF] Zhu | Power-aware scheduling for and/or graphs in real-time systems[END_REF] with AND/OR model applications. Since the target system for both works is shared-memory multiprocessor systems, communication between dependent tasks is not considered unlike our approach.

In [START_REF] Ge | Performance-constrained distributed dvs scheduling for scientific applications on power-aware clusters[END_REF], two voltage scaling algorithms for periodic, sporadic, aperiodic tasks on a dynamic priority single-processor system are proposed. They are more practical compared with many existing DVS algorithms in that a priori information on incoming tasks is not assumed to be available until the tasks are actually released. This assumption does not correspond to our task model as explained.

Rountree et al. in [START_REF] Rountree | Bounding energy consumption in large-scale MPI programs[END_REF] developed a system based on linear programming that exploits slack using DVS (i.e., slack reclamation). Their linear programming system aims to deliver near-optimal schedules that tightly bound optimal solutions. It incorporated allowable time delays, communication slack, and memory pressure into its scheduling. The linear programming system mainly deals with energy reduction for a given pre-generated schedule with a makespan constraint as in most existing algorithms. In our apprche, the makespan is not a constraint but an objective to optimize.

Another two scheduling algorithms for bag-of-tasks applications on clusters are proposed in [START_REF] Kim | Power aware scheduling of bag-oftasks applications with deadline constraints on DVS-enabled clusters[END_REF]. Tasks in a bag-of-tasks application are typically independent and homogeneous, yet run with different input parameters/files. In [START_REF] Kim | Power aware scheduling of bag-oftasks applications with deadline constraints on DVS-enabled clusters[END_REF], deadline constraints are associated with tasks for the purpose of quality control. The two algorithms differ in terms of whether processors in a given computer cluster are time-shared or space-shared. Computer clusters in this paper are composed of homogeneous DVS-enabled processors unlike our approach where processors are heterogeneous.

In [START_REF] Wang | Towards energy aware scheduling for precedence constrained parallel tasks in a cluster with DVFS[END_REF], the authors propose a formulation of energy aware scheduling algorithm and a detailed discussion of slack time computation. This scheduling algorithm also concerns reducing voltages during the communication phases between parallel jobs. In [START_REF] Khan | A cooperative game theoretical technique for joint optimization of energy consumption and response time in computational grids[END_REF], the authors study the energy-aware task allocation problem for assigning a set of tasks onto the machines of a computational grid each equipped with DVS feature. The goal is to optimize the energy consumption and response time in computational grids. Unlike our approach, [START_REF] Khan | A cooperative game theoretical technique for joint optimization of energy consumption and response time in computational grids[END_REF] suppose that tasks are independent and are not subject to precedence constraints.

To the best of our knowledge, none of previous scheduling approaches explicitly addresses the energy issue with a multi-objective approach when tackling the problem of scheduling precedence-constrained parallel applications on HCSs. Therefore, the scheduling algorithms with energy consciousness presented in this section are the most closely related works to our study.

Energy-conscious scheduling heuristic

The consideration of energy consumption in task scheduling adds another layer of complexity to an already intricate problem. Unlike real-time systems, applications in our study are not deadline-constrained. Therefore, the evaluation of the quality of schedules should be measured explicitly considering both makespan and energy consumption. For this reason, energy-conscious scheduling heuristic (ECS) [START_REF] Lee | Minimizing energy consumption for precedence-constrained applications using dynamic voltage scaling[END_REF] is devised with relative superiority (RS) as a novel objective function, which takes into account these two performance criteria.

For a given ready task, its RS value on each processor is computed using the current best combination of processor and VSL (p' and v' are, respectively, the selected processor and its voltage supply level) for that task, and then the processor-from which the maximum RS value is obtained-is selected. Since each scheduling decision that ECS makes tends to be confined into a local optimum, another energy reduction technique (MCER) is incorporated into the energy reduction phase of ECS without sacrificing time complexity. It is an effective technique in lowering energy consumption, although the technique may not help schedules escape from local optima. MCER is makespan conservative in that changes it makes (to the schedule generated in the scheduling phase) are only validated if they do not increase the makespan of the schedule. For each task in a DAG, MCER considers all of the other combinations of task, host and VSL to check whether any of these combinations reduces the energy consumption of the task without increasing the current makespan.

The results clearly demonstrate the superior performance of ECS over DBUS and HEFT. Note that, in many previous studies [START_REF] Kim | Push-pull: Deterministic search-based dag scheduling for heterogeneous cluster systems[END_REF], [START_REF] Lee | A productive duplication-based scheduling algorithm for heterogeneous computing systems[END_REF], HEFT has been proven to perform very competitively, and it has been frequently adopted and extended.

However, ECS returns one solution as a result, and precedence-constrained applications problem is bi-objective in nature. Section 4 presents our new parallel bi-objective approach based on hybridization between a genetic algorithm and ECS. This approach provides a set of solutions to this problem. The experiments presented in Section 5 show that our approach often gives solutions which are better than those found by ECS.

A parallel hybrid approach

This section starts with a brief overview on multi-objective combinatorial optimization and genetic algorithms. Then, our new parallel bi-objective hybrid approach is presented.

Multi-objective Combinatorial Optimization

A multi-objective optimization problem (MOP) consists generally in optimizing a vector of nb obj objective functions F (x) = (f 1 (x), . . . , f nb obj (x)), where x is an d-dimensional decision vector x = (x 1 , . . . , x d ) from some universe called decision space. The space the objective vector belongs to is called the objective space. F can be defined as a cost function from the decision space to the objective space that evaluates the quality of each solution (x 1 , . . . , x d ) by assigning it an objective vector (y 1 , . . . , y nb obj ), called the fitness (Figure 2). While single-objective optimization problems have a unique optimal solution, a MOP may have a set of solutions known as the Pareto optimal set. The image of this set in the objective space is denoted as Pareto front. For minimization problems, the Pareto concepts of MOPs are defined as follows (for maximization problems the definitions are similar):

• Pareto dominance: An objective vector y 1 dominates another objective vector y 2 if no component of y 2 is smaller than the corresponding component of y 1 , and at least one component of y 2 is greater than its correspondent in y 1 i.e.: • Pareto optimality: A solution x of the decision space is Pareto optimal if there is no solution x ′ in the decision space for which F (x ′ ) dominates F (x).

   ∀i ∈ [1..nb obj ], y 1 i ≤ y 2 i ∃j ∈ [1..nb obj ], y 1 j < y 2
• Pareto optimal set: For a MOP, the Pareto optimal set is the set of Pareto optimal solutions.

• Pareto front: For a MOP, the Pareto front is the image of the Pareto optimal set in the objective space.

Graphically, a solution x is Pareto optimal if there is no other solution x ′ such that the point F (x ′ ) is in the dominance cone of F (x). This dominance cone is the box defined by F (x), its projections on the axes and the origin (Figure 3).

Genetic algorithms

Genetic Algorithms (GAs) are meta-heuristics based on the iterative application of stochastic operators on a population of candidate solutions. At each iteration, solutions are selected from the population. The selected solutions are recombined in order to generate new ones. The new solutions replace other solutions selected either randomly or according to a selection strategy.

In the Pareto-oriented multi-objective context [START_REF] Deb | Multi-objective optimization using evolutionary algorithms[END_REF], the structure of the GA remains almost the same as in the mono-objective context. However, some adaptations are required mainly for the evaluation and selection operators.

The selection process is often based on two major mechanisms: elitism and sharing. They allow respectively the convergence of the evolution process to the best Pareto front and to maintain some diversity of the potential solutions. The elitism mechanism makes use of a second population called a Pareto archive that stores the different non-dominated solutions generated through the generations. Such an archive is updated at each generation and used by the selection process. Indeed, the individuals on which the variation operators are applied are selected either from the Pareto archive, from the population or from both of them. The sharing operator maintains the diversity on the basis of the similarity degree of each individual compared to the others. The similarity is often defined as the Euclidean distance in the objective space.

Hybrid approach

In our approach illustrated in Figure 4, a solution (chromosome) is composed of a sequence of N genes. The i th gene of a solution s is denoted s j . Each gene is defined by a task, a processor and a voltage. These three parts of s j are denoted respectively t(s j ), p(s j ) and v(s j ). This means that the task t(s j ) is assigned to the processor p(s j ) with the voltage v(s j ). The new approach we propose is based on ECS [START_REF] Lee | Minimizing energy consumption for precedence-constrained applications using dynamic voltage scaling[END_REF] which is not a population-based heuristic. ECS tries to construct in a greedy way one solution using three components.

• The first component to build the task parts of each gene of the solution.

• The second component to build the processor and voltage parts of these genes.

• And the third component to calculate the fitness of a solution in terms of energy consumption and makespan.

Unlike ECS, our approach provides a set of Pareto solutions. This approach is a hybrid between a multi-objective GA and the second component of ECS.

The role of the GA is to provide good task scheduling. In other words, the GA builds task parts t(s 1 ), t(s 2 ), ..., t(s n ) of a solution s. Therefore, the mutation and crossover operators of the GA affect only the task part of the genes of each solution.

The second component of ECS is called whenever a solution is modified by these two operators. The first role of this component is to correct the task order to take into account the precedence constraints in the task graph. Then the component completes the processor and voltage parts of the genes of the partial solutions provided by these operators. In other words, ECS builds the remaining parts p(s 1 ), p(s 2 ), ..., p(s n ) and v(s 1 ), v(s 2 ), ..., v(s n ) of the partial solutions provided by the mutation and crossover operators of the GA.

The evaluation (fitness) operator of the GA is called once the task, processor and voltage parts of each gene of the solution are known. The role of this operator is to calculate the energy consumption and the makespan of each solution.

The mutation operator is based on the first component of ECS. This first component returns all tasks scheduled according to their b-level values. The principle of our mutation operator is also based on the scheduling of tasks according to their b-level values. The b-level concept is explained in Section 2. It should be noted that one can choose the t-level values instead of those of b-level. First, the operator chooses randomly two integers i and j such that 1 ≤ i < j ≤ n and blevel(t(s i )) < blevel(t(s j )). Then, the operator swaps the two tasks t(s i ) and t(s j ) (see Figure 5).

As illustrated in Figure 6, the crossover operator uses two solutions s1 and s2 to generate two new solutions s ′ 1 and s ′ 2. To generate s ′ 1, the operator:

• considers s1 as the first parent and s2 as the second parent.

• randomly selects two integers i and j such that 1 ≤ i < j ≤ n. • copies in s ′ 1 all tasks of s1 located before i or after j. These tasks are copied according to their positions (s

′ 1 k = s1 k if k < i or k > j).
• copies in a solution s all tasks of s2 that are not yet in s ′ 1. Thus, the new solution s contains (ji + 1) tasks. The first task is at position 1 and the last task at the position (ji + 1).

• and finally, copies all the tasks of s to the positions of s ′ 1 located between i and j (s

′ 1 k = s k-i+1 for all i ≤ k ≤ j).
The solution s ′ 2 is generated with the same method by considering s2 as the first parent and s1 as the second parent. The other elements of the GA in the new approach are conventional. Indeed, our GA randomly generates the initial population. Its selection operator is based on a tournament strategy. The algorithm stops when no new best solution is found after a fixed number of generations.

Insular approach

The island model [START_REF] Cohoon | Punctuated equilibria: A parallel genetic algorithm[END_REF] is inspired by behaviors observed in the ecological niches. In this model, several evolutionary algorithms are deployed to evolve simultaneously various populations of solutions, often called islands. As shown in Figure 7, the GAs of our hybrid approach asynchronously exchange solutions. This exchange aims at delaying the convergence of the evolutionary process and to explore more zones in the solution space. For each island, a migration operator intervenes at the end of each generation. Its role consists to decide the appropriateness of operating a migration, to select the population sender of immigrants or the receiver of emigrants, to choose the emigrating solutions, and to integrate the immigrant ones.

Multi-start approach

Compared to the GA, ECS is more costly in CPU time. The different evaluations of ECS are independent of each other. Therefore, their parallel execution can make the approach faster. The objective of the hybrid approach is to improve the quality of solutions. The island approach also aims to obtain solutions of better quality. The goal of the parallel multi-start approach is to reduce the execution time. As shown in Figure 8, our parallelization is based on the deployment of the approach using the farmer-worker paradigm. The GA processes are farmers and ECS processes are workers.

Experiments and results

This section presents the results obtained from our comparative experimental study. The experiments aim to demonstrate and evaluate the contribution of the hybridization, the insular approach and the multi-start approach respectively compared to ECS, the hybrid approach and the insular approach.

Experimental settings

The performance of our approach was thoroughly evaluated with the Fast Fourier Transformation [START_REF] Cormen | Introduction to algorithms[END_REF] task graph which is a real-world application. A large number of variations were made on this task graph for more comprehensive experiments. Various different characteristics of processors were also applied. Table 4 summarizes the parameters used in our experiments. The new approach is experimented on 9,000 instances distributed equitably according to the number of tasks, the number of processors, the processor heterogeneity and the CCR (1/5 of instances have a number of tasks equal to ∼20, 1/5 of instances have a number of tasks equal to ∼40,..., 1/6 of instances have a number of processors equal to 2, etc.). In other words, 20 instances are used for each combination of parameters.

Our approach has been implemented using ParadisEO [START_REF] Liefooghe | A software framework based on a conceptual unified approach for evolutionary multiobjective optimization: ParadisEO-MOEO[END_REF]. This software platform provides tools for the design of parallel meta-heuristics for multiobjective optimization. [START_REF] Talbi | Metaheuristics: From design to implementation[END_REF] explains how to implement a multi-objective genetic algorithm, an insular approach, and a multi-start using ParadisEO. Table 5 shows the parameters used by ParadisEO for the hybrid, insular and multi-start approaches during our experiments. Experiments have been performed on a grid of three clusters. The first cluster contains 8 Opteron 244 nodes (dual-processor clocked at 1.8 GHz, 2 GB of RAM). The second contains 10 Xeon L5420 nodes (bi-quad-core processors clocked at 2.5 GHz, 16 GB of RAM). The third cluster contains 106 cores AMD opteron 248, 40 cores AMD opteron 252, 104 cores AMD opteron 285, et 368 cores Intel Xeon E5440 QC. A total of 714 cores are used. The first two clusters are located at the University of Mons in Belgium, while the third cluster is at Université de Lille1 in France.

Hybrid approach

The hybrid approach is experimented on all instances of Table 4. Each instance is solved twice. The first resolution is done with ECS, and the : An example of the obtained results with the hybrid approach and ECS for the same instance second resolution with the new approach. These experiments are launched by a script on one of the cores of our grid according to their availability. Figure 9 gives an example of a Pareto front obtained with the hybrid approach and the solution obtained by ECS for the same instance which is the tenth instance generated with the number of tasks equal to ∼20, number of processors equal to 02, processor heterogeneity equal to 100, and CCR equal to 0.1. Experiments show that ECS finds the solution of an instance after ∼ 1 second on average, while our hybrid approach requires about ∼ 25 minutes on average to find the Pareto solutions of an instance. As previously mentioned, ECS is a heuristic that builds one and only one solution using a greedy strategy, and our hybrid approach is based on the hybridization between a GA and ECS. Therefore, the hybrid approach handles a population of solutions that evolves over several generations, and the second component of ECS is called and used during the construction of each solution. So it was expected that the hybrid approach uses more computing power than ECS. Our goal is not to have an approach faster than ECS but an approach that gives Pareto solutions which improve the solution of ECS. Our approach is useful, for example, for large scientific applications requiring high computing power, and for small applications which are executed many times.

Table 6 compares the Pareto solutions of the hybrid approach with the solution of ECS. The comparison is made according to the number of tasks, the number of processors, the processor heterogeneity and the CCR. The third column shows the average number of obtained Pareto solutions. The last column gives the percentage of Pareto solutions that improves the ECS solution on the two objectives simultaneously. As indicated in the last line of the table, the hybrid approach provided 19.77 solutions on average, and 83.04% of the Pareto solutions found by this hybrid approach improve the ECS solution on the two objectives simultaneously. In addition, Table 6 shows that the more tasks there are, the more Pareto solutions are found, and the more the percentage of Pareto solutions dominating the ECS solution increases. To determine the contribution of the new approach, in terms of the values of makespan and energy consumption, we compare the solution provided by ECS to only one solution of the Pareto set provided by the new approach. The solution chosen in the Pareto set is used only to compare the new approach with ECS. Nevertheless the decision maker, using the new approach, will have a set of Pareto solutions instead of one solution.

For each instance,

• a first resolution is done using ECS to provide one solution s.

• a second resolution is done using the new approach to obtain a set E of Pareto solutions.

• only one Pareto solution s ′ is selected from the set E. This solution is the closest to s in the sense of Euclidean distance.

• finally, a comparison will be done between the solutions s and s ′ .

Figure 10, Figure 11, Figure 12 and Table 7 allow to compare in a detailed way the two approaches. They respectively show the improvement brought by the new approach according to the number of tasks, the number of processors, the CCR, and the processor heterogeneity. Experiments show that our approach improves on average the results obtained by ECS. Indeed, as shown in Table 7, the energy consumption is reduced by 47.49% and the makespan reduced by of 12.05%. In addition, Figure 11 shows clearly that the more processors there are, the more the new approach improves the results of ECS.

Insular approach

The objective of the following experiments is to show that our island approach improves the quality of the solutions provided by the hybrid approach. This insular approach is useful when solving large instances. Therefore, the 4. The instances used are those with the number of tasks is 120, the number of processors is 64, the value of CCR is 10, and the heterogeneity of processors is 200 (20 instances). Each instance is solved using 1, 5, 10, 30 or 50 islands. An insular approach with 1 island is equivalent to the hybrid approach.

Figure 13 illustrates the S-metric average values obtained with different numbers of islands. These values are normalized with the average value obtained by the experiments using 1 island. The S-metric measures the hyper-volume defined by a reference point and a Pareto front. It allows to evaluate the quality of a Pareto front provided by an algorithm.

Experiments show that whatever the number of used islands the insular approach improves the Pareto front obtained with the hybrid approach. As shown in Figure 13, the use of 50 islands, instead of 1 island (i.e. the hybrid approach), improves the S-metric of the obtained Pareto front by 26%. In Figure 13, the more the number of islands is used, the better the results will be.

Multi-start approach

This section presents the experiments done to assess the quality of our multi-start approach. The parameters of the instances used in our experiments are: the CCR is 0.1, 0.5, 1.0, 5.0 or 10.0, the number of processors is to run the GA). In our case, an experiment can not have a speedup greater than 21. Table 8, Figure 14, Figure 15 show respectively the evolution of the speedup according to the processor heterogeneity, the CCR, and the number of processors. The average speedup obtained is 13.06.

As shown in Figure 15, the speedup increases proportionally to the number of processors on which the precedence-constrained parallel application is run. Table 8 and Figure 15 that show the CCR and the heterogeneity of processors do not impact significantly the quality of the acceleration of our approach.

Conclusions

In this paper, we have investigated the precedence-constrained parallel applications particularly on high-performance computing systems like cloud computing infrastructures. Precedence-constrained parallel applications are designed mostly with the sole goal of minimizing completion time without paying much attention to energy consumption. We presented a new parallel bi-objective hybrid genetic algorithm to solve this problem. The algorithm minimizes energy consumption and makespan. The energy saving of our approach exploits the dynamic voltage scaling (DVS) technique-a recent advance in processor design.

Our new approach has been evaluated with the Fast Fourier Transformation task graph which is a real-world application. Experiments show that our bi-objective meta-heuristic improves on average the results obtained in the literature (see [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF], [START_REF] Bozdag | A task duplication based bottom-up scheduling algorithm for heterogeneous environments[END_REF] and [START_REF] Garey | Computers and intractability: A guide to the theory of NP-completeness[END_REF]) particularly in energy saving. Indeed, the energy consumption is reduced by 47.5% and the completion time by 12%. The experiments of the insular approach also show that the more the number of islands is used, the better the results will be. The use of 50 islands, instead of 1 island (i.e. the hybrid approach), improves the S-metric of the obtained Pareto front by 26%. Furthermore, the multi-start approach is on average 13 times faster than the island approach using 21 cores.

However, we observed that the hybrid approach consumes more resources than ECS, and the insular approach consumes more resources than the hybrid approach. In the insular approach, experiments show that the more the number of islands is used, the more the resources are needed. A resource can be a processor, a network bandwidth, etc. The energy consumed by an approach increases when the used resources increase. We think that the multi-start approach does not increase significantly the energy consumed by the insular approach.

Therefore, one of the main perspectives of the work presented in this paper is to determine the solving approach to choose among ECS, the hybrid approach, and the insular approach, according to the precedence-constrained parallel application at hand. If the insular approach is chosen, the major issue is to determine the best number of islands to be used. This future work aims to minimize the total amount of consumed energy by the chosen solving approach and by the precedence-constrained parallel application to be solved. It is clear, for example, that the insular approach is interesting for the large and resource consuming precedence-constrained parallel applications and the applications intended to be executed several times.

Appendix Algorithm 1

The main parameters of our approaches task(solution buffer m )←task(parent2 k )

Figure 1 :

 1 Figure 1: A simple task graph

Figure 2 :

 2 Figure 2: Illustration of a MOP

1 Figure 3 :

 13 Figure 3: Example of non-dominated solutions

Figure 4 :

 4 Figure 4: Our hybrid GA (GA and ECS)

Figure 5 :

 5 Figure 5: The mutation operator

Figure 6 :

 6 Figure 6: The crossover operator

Figure 7 :

 7 Figure 7: The cooperative island approach

Figure 8 :

 8 Figure 8: Illustration of the multi-start approach

Figure 9

 9 Figure 9: An example of the obtained results with the hybrid approach and ECS for the same instance

22

 22

Figure 10 :

 10 Figure 10: Improvement according to the number of tasks

Figure 11 :

 11 Figure 11: Improvement according to the number of processors

Figure 12 :

 12 Figure 12: Improvement according to the CCR

Figure 13 :

 13 Figure 13: S-metric value according to the number of islands

Figure 14 :

 14 Figure 14: Speedup according to the CCR

Figure 15 :

 15 Figure 15: Speedup according to the number of processors

Algorithm 2

 2 The main algorithm of our approaches hybrid insular multistart approches() 1: initialize(population,POPULATION SIZE) 2: evaluate(population) 3: GENERATION←1 4: while GENERATION ≤ GENERATION MAXIMUM do

1 : 1 : 29 Algorithm 5

 11295 for all solution ∈ population do 2: launch parallel(ECS component2,solution) 3: end for 4: for all solution ∈ population do 5: cost(solution)←read cost(solution) 6: end for evaluate sequential(population) 1: for all solution ∈ population do 2: cost(solution)←ECS component2(solution) 3: end for Algorithm 4 Mutation operator mutation(children) 1: for all solution ∈ children do 2:if random([0,1]) ≤ MUTATION RATE then (i,j)=random(1≤i<j≤N ∧ check level(solution,i,j)) 2: swap(task(solution i ),task(solution j )) check level(solution,i,j) 1: return b-level(task(solution i ))<b-level(task(solution j )) Crossover operator crossover(parents) 1: children← ∅ 2: for all i ∈ (1...POPULATION SIZE) do 3:if random([0,1]) ≤ CROSSOVER RATE then end for 10: return children crossover2(parent1,parent2)1: child1←crossover1(parent1,parent2) 2: child2←crossover1(parent2,parent1) 3: return (child1,child2) crossover1(parent1,parent2)1: (i,j)=random(1≤i<j≤N) 2: for all k ∈ (1,...,i-1,j+1,...,N) do 3: task(child k )←task(parent1 k ) 4: end for 5: m←0 6: for all [k ∈ (1,..,N)] ∧ [task(parent2 k ) / ∈ tasks(child)] do 7:

Table 1 :

 1 Voltage-relative speed pairs

		Pair 1	Pair 2	Pair 3
	Level Voltage Relative Voltage Relative Voltage Relative
		(v k )	speed	(v k )	speed	(v k )	speed
			(%)		(%)		(%)
	0	1.5	100	2.2	100	1.75	100
	1	1.4	90	1.9	85	1.4	80
	2	1.3	80	1.6	65	1.2	60
	3	1.2	70	1.3	50	0.9	40
	4	1.1	60	1.0	35
	5	1.0	50
	6	0.9	40

Table 2 :

 2 Computation cost with VSL 0

	Task p 0 p 1 p 2
	0	11 13 9
	1	10 15 11
	2	9 12 14
	3	11 16 10
	4	15 11 19
	5	12 9 5
	6	10 14 13
	7	11 15 10

Table 3 :

 3 Task Priorities

	Task b-level t-level
	0	101.33 0.00
	1	66.67 22.00
	2	63.33 28.00
	3	73.00 25.00
	4	79.33 22.00
	5	41.67 56.33
	6	37.33 64.00
	7	12.00 89.33

Table 4 :

 4 Experimental parameters

	Parameter	Value
	The number of tasks	∼20 ∼40 ∼60 ∼80 ∼120
	The number of processors	02 04 08 16 32 64
	Processor heterogeneity	100 200 random
	CCR	0.1 0.2 1.0 5.0 10.0

Table 5 :

 5 The parameters used by ParadisEO for each approach

	Parameters	Hybrid	Insular	Multi-start
	Population size	20	20	20
	Number of generations	1000	100	100
	Crossover rate	1	1	1
	Mutation rate	0.35	0.35	0.35
	Migration topology		Ring
	Migration rate		Every 20 generations
	Number of migrants		5

Table 6 :

 6 Comparison of Pareto hybrid approach solutions and ECS solution

			Average	Pareto solutions
			number of	dominating ECS
			Pareto solutions	solution (%)
		∼20	14.78	78.24
	Number	∼40	19.57	80.70
	of tasks	∼60	21.36	83.62
		∼80	21.45	83.12
		∼120	21.67	89.51
		02	18.51	73.21
		04	19.42	71.01
	Number	08	22.17	75.83
	of processors	16	23.32	86.12
		32	19.98	94.73
		64	15.18	97.36
		100	20.12	80.82
	Heterogeneity	200	21.67	74.47
		random	17.52	93.83
		0.1	19.60	88.44
		0.2	19.47	88.83
	CCR	1.0	17.53	89.09
		5.0	19.26	78.80
	Average		19.77	83.04

Table 7 :

 7 Improvement according to the processor heterogeneity

	Processor heterogeneity Energy (%) Makespan (%)
	100	44.74	9.10
	200	43.49	7.07
	random	54.26	19.99
	Average	47.49	12.05

Table 8 :

 8 Speedup according to the processor heterogeneity

	Processor heterogeneity Speedup
	100	13.74
	200	12.73
	random	12.73
	Average speedup	13.06

  end for 10: for all k ∈ (i,...,j) do DESTINATION←destination(topology)9: send(DESTINATION,migrants) 10: SOURCE←source(topology) 11: migrants←receive(SOURCE) 12: insert(population,migrants)

	Algorithm 6 Migration operator
	migrate(population,GENERATION)
	1: if NOT INSULAR then
	2:	stop
	3: end if
	4: if GENERATION mod(MIGRATION RATE) =0 then
	5:	stop
	6: end if
	7: migrants←select(population,MIGRANTS SIZE)
	8:
	8:	m←m+1
	9: 11:	task(child k )←task(solution buffer k-i+1 )
	12: end for
	13: return child

Acknowledgments

Professor Albert Zomaya's work is support by an Australian Research Grant DP1097110. Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr). Clusters of University of Mons were also used. We would like to thank the technical staffs of the Grid'5000 and the clusters of University of Mons for making their clusters accessible and fully operational.


F.-M Birleanu

C Ioana

C Gervaise

J Chanussot

A Serbanescu

G Serban

ON THE RECURRENCE PLOT ANALYSIS METHOD BEHAVIOUR UNDER SCALING TRANSFORM

Keywords: Transient signals, recurrence plots, amplitude and time scaling

In the last decade, the applications of the recurrence plot analysis method make it a valuable alternative to the time-frequency and time-scale tools. As it was initially developed for the study of dynamical systems, and was later used in nonlinear time series analysis, the question of using it as a signal processing tool has not been put into discussion yet. In this field the projective techniques are largely used, with good results. Nevertheless, they also have some limitations -especially regarding transient signal processing. But this kind of signals are ubiquitous in real world. In addition, propagation through various media as well as on multiple paths lead to delayed, attenuated and dilated versions of the original transients. In this paper we study the behaviour of the recurrence plot analysis method in the context of analyzing some finite duration signals being subject to rescalings of the amplitude and time axes. This study is a starting point in employing the analysis of recurrences in investigations of a large class of real world signals.

INTRODUCTION

Recurrence plot analysis [1] has its roots in a 1987 paper by Eckman et al [2]. They introduced the recurrence plot as a "new graphical tool for measuring the time constancy of dynamical systems". Although its initial purpose was to obtain dynamical parameters from time series, it was subsequently applied in various areas of the field of nonlinear time series analysis. According to [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF], the methods based on recurrence plots have been successfully applied in "physiology, neuroscience and genomics, ecology, physics, chemistry, earth science and astrophysics, engineering and economy". The bibliography gathered on the recurrence plot website (http://www.recurrence-plot.tk) shows significant increasing of interest for this tool in the last ten years, with a peak in the last three years. However, very few applications of recurrence plots were reported in signal processing -mostly in speech processing and signal detection (e.g. [START_REF] Schinkel | Selection of recurrence threshold for signal detection[END_REF][START_REF] Michalowicz | Signal detection based on recurrence matrix statistics[END_REF][START_REF] Rohde | Stochastic analysis of recurrence plots with applications to the detection of deterministic signals[END_REF][START_REF] Marwan | Comment on "stochastic analysis of recurrence plots with applications to the detection of deterministic signals" by rohde et al[END_REF]). The integration of this technique in the signal processing field seems natural, as there is a strong similarity between recurrence and frequency. This represents the main objective of our paper. Hence, we start by studying the effects of scaling the amplitude and the time axis of the signal on the resulting recurrence plot. Scaling phenomena (like attenuations and dilatations) are common in real world signals. The issue of dilatations has been solved by the wavelet transform, but the choice of the proper mother wavelet is not always possible for all signals. On the other hand, recurrence plot analysis does not need a reference signal, as it only compares the analyzed signal to itself.

The paper is organized as follows. In the following, we introduce first the principle of recurrence plot analysis and we continue by presenting the effects of scaling the signal amplitude. Then, we present the effects of time scaling and we continue with a short discussion regarding the potential of the obtained results. We close by pointing out the conclusions.

RECURRENCE PLOT ANALYSIS

Figure 1 summarizes the generic steps of the recurrence plot analysis method. They are briefly discussed in the following paragraphs.

(1) -Representation of the signal as a trajectory in a multidimensional space. This is performed with the aid of the method of delays. If m is the dimension of the representation space and τ is the time delay, then the phase space trajectory corresponding to the signal s(t) is:

r(t) = m X k=1 s (t + (k -1)τ ) • e k , (1)
where e k are the versors of the phase space axes. Hence, the time evolution of the trajectory will be the same on all the m axes, except for time delayings by multiples of τ .

(2) -Computation of the recurrence plot. A recurrence between points at times i and j on the trajectory is defined as:

R(i, j) = Θ (ε(i) -D ( r(i), r(j))) , (2)
where Θ denotes the Heaviside step function, and D ( r(i), r(j)) stands for the distance between points i and j on the trajectory. In fact, Equation (2) shows that a recurrence is identified between points i and j on the trajectory if the distance between them is smaller than a threshold ε(i), known as the recurrence radius. Most often, a constant ε is used, and the distance D is computed using the Euclidean metric. (We noticed that in most cases it is sufficient to choose ε as the mean distance between successive points of the trajectory.) However, it might sometimes be useful to work directly with the distance plot (i.e the graphical representation of D).

(3) -Quantification of the recurrence plot. It consists in performing different computations on the recurrence plot, in order to obtain some measures that are able to offer some insight into the analyzed signal. The most common such measure is the recurrence (1) (2)

(3) (3) quantification of the recurrence plot. In this example, signal s(t) is obtained from a chaotical Lorenz system. The trajectory r(t) is obtained using a dimension m = 3 and a delay τ = 3. It can be noticed that r(t) has strong similarities with the well-known Lorenz "butterfly" attractor. The transitions of the system from one unstable periodic orbit to the other can be identified using the significant maxima of the locally-computed (by using a window of size 10) recurrence rate, RR(t).

rate, RR. It is computed as the ratio between the total area occupied by black dots (recurrences) in the recurrence plot and the total area of the recurrence plot. The recurrence rate, as well as other quantification measures, can be computed either globally, for the whole recurrence plot, or locally, on time windows that are being shifted along the main diagonal of the recurrence plot. In this way, a timevarying measure, RR(t), will be obtained. Various recommendations can be found in literature regarding the choice for the parameters of the recurrence plot analysis method (i.e. m, τ , ε, D) [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF]. Only the most commomly used recurrencebased analysis method was presented in this section.

AMPLITUDE SCALING

As the phase space representation of the signal is a linear operation, scaling the signal amplitude will result in an appropriate scaling of the trajectory. If for signal s(t) we have trajectory r(t), then for signal s ′ (t) = α • s(t) we have:

r ′ (t) = m X k=1 α • s (t + (k -1)τ ) • e k = α • r(t). (3)
Let us see whether this scaling of the trajectory has any effect on the recurrence plot. If the recurrence plot corresponding to signal s(t) is R(i, j), then the recurrence plot corresponding to signal s ′ (t) will be:

R ′ (i, j) = Θ `ε′ (i) -D (α • r(i), α • r(j)) ´. (4)
The distance D is usually computed using a certain metric, that is:

D ( r(i), r(j)) = r(i) -r(j) . (5)
In this case, in order for R ′ (i, j) to be identical to R(i, j), ε ′ (i) should be chosen such that:

ε ′ (i) = α • ε(i). ( 6
) (a) (b)
Fig. 2. The shape of the neighbourhood (in a tridimensional phase space) is a sphere for the Euclidean metric (a) and a cone for the angular distance (b).

However, this scaling of the recurrence radius is no more needed when using an angular distance, defined as:

D ( r(i), r(j)) = arccos " r(i) • r(j) r(i) • r(j) « , (7)
for non-zero vectors (otherwise vectors can be considered paralell, i.e. the angular distance between them is 0). This definition for D leads to a recurrence plot that is invariant to amplitude scalings of the analyzed signal. We note that this angular distance generates conical neighbourhoods (instead of the spherical ones generated by the Euclidean metric). Figure 2 illustrates this observation. We also note that by using such a conical neighbourhood the concept of recurrence loses its initial meaning (that is returning close to a previously visited point) -therefore, we are actually working with a generalized recurrence.

Fig. 3. Time-scaled one-period sinusoidal signal and the corresponding trajectories. From top to bottom, β is: 1, 1.73, 0.61, and τ is: 250, 250/1.73, 250/0.61. The shape of the trajectory remains unchanged.

TIME SCALING

If s(t) is the signal, the time-scaled version of it is s ′ (t) = s(β • t).

The corresponding trajectory will be:

r ′ (t) = m X k=1 s `βt + (k -1)βτ ′ ´• e k . (8)
We notice that scaling the time axis of the signal leads to an appropriate scaling of the time delay (the τ parameter). For the trajectory to remain unchanged, τ ′ should be chosen such that:

τ ′ = τ β . (9)
Figure 3 shows that the trajectories indeed remain unchanged when τ is properly scaled. However, Equation [START_REF] Marwan | Encounters with neighbours. Current developments of concepts based on recurrence plots and their applications[END_REF] shows that r ′ (t) is not identical to r(t). Figure 3 does not reflect the time evolution, and therefore r ′ (t) and r(t) seem identical. In fact, trajectory r(t) is "drawn" β times faster (or slower, depending on whether β is greater than or less than 1) than trajectory r ′ (t). This leads to the contraction/dilatation of the structures in the recurrence plot, as shown in Figure 4.

DISCUSSION

We have shown in section 3 that we can obtain a recurrence plot that is invariant to amplitude scaling by using the angular distance. As for the time scaling, we showed in the previous section that by appropriately scaling the time delay the shape of the trajectory remains unchanged (Fgure 3), except it is covered with a different speed, which results in a stretching of the recurrence plot pattern (Figure 4). Failure to meet the constraint in Equation ( 9), though, leads to As τ is scaled appropriately, the pattern remains unchanged, except for the fact that it is stretched.

significant changes in the trajectory -redundance (when β < 1) or irrelevance (when β > 1) [START_REF] Marwan | Encounters with neighbours. Current developments of concepts based on recurrence plots and their applications[END_REF]. It follows that the recurrence patterns will also be modified. However, when dealing with dilated signals, i.e. when β < 1, the trajectory does not change its main characteristics, except that it is "attracted" to the first diagonal of the coordinate system. (In this case the changes in the recurrence pattern are mostly due to the use of recurrence balls that have the same dimensions as those used for the non-dilated signal.) Figure 5 shows an example of an angular distance plot computed for a synthetic signal (and Figure 6 repeats the example, after adding some noise). The signal is composed of four transients coming from two different sources -the first two of them come from direct propagation between the sources and the receiver, and the last two come from reflexions (i.e. they are attenuated and dilated). As the figure shows, the angular distance plot of this signal allows a good time localization of the four transients, regardless of their scale differences. More than that, it allows a visual characterization of the transients. It can be easily observed (by analyzing the blocks on the main diagonal, that correspond to the four transients) that we are dealing with two types of signals (i.e. two different sources). It can also be easily observed that the first transient is of the same type as the third, and that the second is of the same type as the fourth. We can make this statement either as a result of a comparative visual analysis of their corresponding blocks in the angular distance plot, or as a result of the fact that their cross-term blocks are composed of diagonally oriented structures (e.g. the third rectangular block on the first line of D(i, j) in Figure 5 is the cross-term block between the first and the third transients; all the black structures in it are oriented along the first diagonal of the block).This final observation may offer, together with the observations made in section 4, a good starting point for developing a robust automatic method for comparing transient signals.

CONCLUSION

We presented in this paper the essential aspects of a nonlinear time series analysis method called recurrence plot analysis. As this tool is working its way into signal processing, we proposed to study its properties in the context of modifying the amplitude scale as well as the time scale of the signal. We showed that an amplitude scaling does not lead to a different recurrence plot as long as the recurrence radius is scaled appropriately. Besides that, we showed an alternative method to compute the recurrence plot, by using the angular distance between vectors, and we showed that the recurrence plot obtained this way is invariant to signal amplitude scaling.We also showed that a time scaling of the signal does not lead to a modified trajectory if the time delay is appropriately scaled. The recurrence plot remains also unchanged, except for the fact that it will compact/expand (without otherwise affecting the characteristics of the patterns it contains). This theoretical study may serve as a basis for involving recurrence plot analysis in solving problems in the analysis of transient signals, that are encountered in many real world applications, where attenuations and dilatations are very common. Further, we are planning to use the observations made in this paper for developping a scale-independent method for the characterization of transient signals.

[1] N. Marwan, "A historical review of recurrence plots," The European Physical Journal Special Topics, vol. 164, pp. [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF][START_REF] Schinkel | Selection of recurrence threshold for signal detection[END_REF][START_REF] Michalowicz | Signal detection based on recurrence matrix statistics[END_REF][START_REF] Rohde | Stochastic analysis of recurrence plots with applications to the detection of deterministic signals[END_REF][START_REF] Marwan | Comment on "stochastic analysis of recurrence plots with applications to the detection of deterministic signals" by rohde et al[END_REF][START_REF] Marwan | Encounters with neighbours. Current developments of concepts based on recurrence plots and their applications[END_REF][9][10][11][12]2008.

[2] J.-P. Eckmann, S. Oliffson Kamhorst, and D. Ruelle, "Recur- rence plots of dynamical systems," Europhysics Letters, vol. 4, no. 9, pp. 973-977, November 1987.

Fig. 1 .

 1 Fig. 1. The steps of recurrence plot analysis: (1) phase space representation of the signal; (2) computation of the recurrence plot;(3) quantification of the recurrence plot. In this example, signal s(t) is obtained from a chaotical Lorenz system. The trajectory r(t) is obtained using a dimension m = 3 and a delay τ = 3. It can be noticed that r(t) has strong similarities with the well-known Lorenz "butterfly" attractor. The transitions of the system from one unstable periodic orbit to the other can be identified using the significant maxima of the locally-computed (by using a window of size 10) recurrence rate, RR(t).

Fig. 4 .

 4 Fig. 4. The effect of dilatation on the recurrence plot, for a sinusoidal signal. (left) The recurrence plot of signal s(t) is obtained using m = 2 and τ = 50. (right) The recurrence plot of signal s(0.61 • t) is obtained using m = 2 and τ ′ = 82(= τ /0.61). As τ is scaled appropriately, the pattern remains unchanged, except for the fact that it is stretched.

Fig. 5 .

 5 Fig. 5. Angular distance plot for a synthetic signal composed of four transients coming from two different sources. (The signal-to-noise ratio is 60 dB.) It is computed using m = 3 and τ = 15. The four transients can be easily localized in time and characterized as well by a visual analysis of this image.

Fig. 6 .

 6 Fig.6. Angular distance plot for the same signal in Figure5, with a signal-to-noise ratio of 5 dB. It is computed using the same parameters (m = 3 and τ = 15). Noise has a strong impact on the resulted image, but time localization and characterization of the four transients can still be done to some extent.


F.-M Birleanu

-C Ioana

- A Serbanescu

G Serban

-E Sofron

Estimation of thermo-hydrodynamic parameters in energy production systems using non-stationary signal processing

Keywords: Non-intrusive measurement, ultrasounds, temperature and flow rate estimation, polynomial phase modeling. I

In order to control water flow in power generation systems, obtaining precise information about the parameters of the flow is important. The construction of such systems, as well as the operating conditions, impose the use of non-intrusive techniques. The use of ultrasounds allows solving this problem in a simple and elegant manner. The technique consists in sending a signal to be affected by the flow without sensibly interacting with it, and then extracting the desired parameters of the flow from the received signal. The paper first presents how the acoustic wave is altered by the propagation in the liquid environment, and then it studies the means to analyze the received signal in order to obtain from it information about two important parameters of the flow -the temperature of the water and the flow rate. Polynomial phase modeling is used, for obtaining a good estimation of the propagation times.

INTRODUCTION

The generation of electrical energy involves water flows through pipes. The parameters of these flows may provide information about the state of the system or of the pipe. In addition to that, the well-functioning of the system implies mantaining these parameters in between certain limits. Therefore, these flows must be controlled, and in order to close the loop of the control system, a method to obtain information about the flow is needed. This method is constrained by the fact that it must work while the system is on, and without affecting its functioning. A non-intrusive method is, therefore, needed. Ultrasounds promise to offer an elegant and efficient solution. In brief, the method involves the following: mounting the transducers on the outer walls of the pipe, transmitting a signal to be affected by the flow, and analyzing the received signal in order to obtain the flow parameters.

In this paper we will deal with the estimation of two very important parameters of the water inside the pipes -its temperature and its flow rate. We first show how the ultrasonic wave is modified by the interaction with the environment, in the contexts of determining each of these two parameters. We then present how the received signal is analyzed in the purpose of obtaining the information that allows estimating the desired thermo-hydrodynamic parameters.

The rest of the paper is organized as follows. Section 2 briefly presents the ultrasonic waves and the way they are affected by the temperature of the water through which they propagate, as well as by its flow rate. Section 3 presents the methods that were used to extract the desired information from the received signals. In section 4 we present the results we obtained using the described materials and methods. We conclude the paper with section 5, in which we discuss the results.

II. ULTRASOUND PROPAGATION THROUGH WATER

The ultrasonic wave [START_REF] Papadakis | Ultrasonic instruments and devices[END_REF] is a pressure wave that propagates through an elastic medium with a frequency that is above 20 KHz. That means that is consists in the propagation of a mechanical energy in a medium. The propagation medium of the ultrasonic wave is subject to a succession of over-pressures an under-pressures, and its constitutive particles are then animated by a back and forth movement on the propagating direction of the wave, as illustrated in Figure 1. As any wave, the ultrasonic wave will interact with the propagation medium, being subject to attenuations, reflexions, refractions and diffractions.

We must note that before propagating through water the wave must first pass through the steel walls of the pipe. Ultrasound wave propagation through such solid, non-elastic media is possible, but is performed as transversal waves, with a much higher speed.

II.1 Measuring temperature

Sound propagation speed in pure water depends on some properties of the water, such as its local density [START_REF] Raichel | The Science and Applications of Acoustics[END_REF], [START_REF] Everest | The Master Handbook of Acoustics[END_REF]. This speed can be found starting from the state equation of pure water. However, there are many simplified versions. The evolution law of the sound speed can be reduced to a dependency on the temperature of the fluid and, less importantly, to its pressure. The following equation presents the evolution of the sound speed in pure water with respect to the temperature, and it is obtained by polynomial interpolation of the experimental results: ! c = 1.4024 "10 3 + 5.0384 " T # 5.8114 "10 #2 " T 2 + 3.3456 "10 #4 " T 3 #1.4815 "10 #6 T 4 + 3.1608 "10 #9 " T 5

.

(1)

Graphically, this evolution looks like in Figure 2. This strong dependence of the speed of sound with temperature will be used to determine the latter. There is, however, a problem. As seen in Figure 2, for a certain range of values for ! c , there are two possible values for ! T . Therefore, is seems impossible to determine ! T knowing ! c . One possible solution to this problem would be estimating the level of noise in the received signal, as it is related to the temperature -a high level of noise would indicate a greater temperature.

II.2 Measuring flow rate

If an ultrasonic impulse is sent to a particle that is at rest, the echo produced by this impulse will be of the same shape. If the particle gets off the emitter, the ultrasonic wave will need a longer time to cover the way to the particle and back. This variation of the length of the route, caused by the variation of the distance to be covered, is called propagation delay. The echos of a particle are always the same when the particle is at rest, that is, when there is no propagation delay. The echos have the same relative phase, which means that there is no phase difference. In case of superposition of a second echo on a particle that gets off the receiver, this echo will be delayed with respect to the first one.

A propagation delay corresponds to a modification of distance, so if we measure the delay and knowing the speed of sound, we can evaluate how far the particle moved, and as we know the delay between the pulses, we can compute the speed of the particle [START_REF] Lynnworth | Ultrasonic Measurements for Process Control: Theory, Techniques, Applications[END_REF]. This measurement principle is illustrated in Figure 3. The flow rate measurement principle using the difference of propagation times consists in performing transmissionsreceptions in the flow direction, as well as in the reverse direction. In the first case, the speed of fluid particles will affect the covered time, ! T AB , as follows:

! T AB = L c + v cos" + 2T paroi + 2T capteur , (2)
where ! c is the speed of sound, ! v is the average speed of the fluid, ! " is the angle between the propagating direction and the horizontal, ! T paroi is the propagation time through the walls and ! T capteur is the transit time inside the transducers. We may note that the speed of the fluid acts as an additional apparent speed that is added when the emission is performed in the direction of the flow. In the other case, when the emission is performed in the reverse direction, this apparent speed diminishes the global speed of sound. All these lead to a propagation time like:

! T BA = L c " v cos# + 2T paroi + 2T capteur . (3)
Using these two measures of the arrival time, the speed of the fluid can be expressed as:

! v = L 2 2D "T t AB # t BA $ % & ' ( ) , with ! "T = T BA # T AB , (4)
where ! t AB = T AB " 2T paroi " 2T capteur and ! t BA = T BA " 2T paroi " 2T capteur . The volumic flow rate is then computed starting from the average speed, using:

! Q = "D int 2 4 # L 2 2D # $T t AB # t BA % & ' ( ) * # 1 K H , (5)
where ! D int is the inner diameter of the cylinder and ! K H is the hydraulic coefficient, computed starting from the Reynolds number [START_REF] Lesieur | Turbulence in Fluids[END_REF], Re, using:

! K H = 2n + 1 n , with ! 1 n = 0.25 " 0.023# log Re ( ) . (6)
What has been discussed so far in this section shows that estimating the propagation times gives us access to the flow rate of the fluid, as well as to its temperature. We will show in the next section a signal analysis method that is focused on estimating as accurate as possible these arrival times.

II.3 Experimental setup

The block diagram of the experimental setup used for measuring the temperature is presented in Figure 4.

!"#$"%&'(%" $%)*" )% +),%&-". /"-"%&')% 0123..)12)$" 4+563-'"%,&2"6! ')6&-&.71316)-689 We used hot tap water, and we let it get cold progresively. The non-intrusive temperature measurement procedure was repeated every ten minutes. The temperature was also measured using a thermometer that was introduced inside the water, at the level of the transducers.

The flow rate was measured using an experimental setup whose block diagram in presented in Figure 5. Several positions for the transducers have been tried, in order to maximize the received energy, using the available equipment ( -the sampling frequency of the oscilloscope was 10 MHz).

III. THE ANALYSIS OF THE RECEIVED SIGNALS

A brute visual analysis of the received signal shows that before the expected signal we obtain a signal that corresponds to the propagation (with a speed that is much greater than the speed in water) through the steel walls of the pipe. The fact that this separation and detection of the signals can be done visually suggests that we could analyze them using methods that are based on the analysis of the temporal shape of the received pulses. We must note, however, that noise and distorsions due to the propagation pose problems to this class of methods.

A possible solution consists in considering the envelope of the received pulses as a suite of phase modulations and trying to estimate them by a polynomial phase modeling. This concept, introduced in [START_REF] Ioana | Time-Frequency Analysis using Warped-Based High-Order Phase Modeling[END_REF], has the advantage of providing a representation space that is more robust to perturbations of any kind: noise, propagation, multiple routes, etc. More precisely, the envelope of the received pulses is characterized by the polynomial coefficients whose values stay stable in the presence of perturbations.

Mathematically, a warping operator is defined as a unitary transform U on l 2 (!), whose effect is x-axis deforming : ,

where ! w is a smooth, one-to-one function, comprising a large subclass of unitary transformations. Generally, these functions, called warping functions, are chosen to ensure the "linearization" of signal time-frequency behavior. To reduce the error propagation effect, the warping operator is defined as as :

,

where is the k th order estimated polynomial coefficient.

Polynomial modeling is done using higher order ambiguity functions [START_REF] Ioana | Time-Frequency Analysis using Warped-Based High-Order Phase Modeling[END_REF]. In our context, this modeling starts with evaluating the ambiguity function of order 6 that allows, through an estimation of maxima, the estimation of the polynomial coefficient of order 6. Its influence on the analyzed signal is removed using a warping procedure, thus avoiding the propagation of the estimation errors that might appear, to lower orders. This procedure, that is called WarpCom, is illustrated in Figure 6 and is followed until the estimation of the first order coefficient is done. For the pulse train in Figure 7, the polynomial modeling of the emitted and, respectively, of the received signal are illustrated in the same figure. We may note the quasiinvariance we obtained thanks to the polynomial coefficients. In order to reduce the effect of modeling errors, we may proceed to an averaging on more pulses of the results of the polynomial modeling.

IV. RESULTS

The method we presented offers good results for estimating the speed of sound in terms of temperature. The results obtained by estimation are illustrated comparatively together with the thoretical ones in Figure 8. As the desired parameter is in fact the temperature, it can be found starting from the speed of sound and using the curve in Figure 2. As mentioned before, for a certain range of speeds there is no bi-univocal correspondence between speed and temperature -thus an estimation of the level of noise (that increses with temperature) might discriminate between the two cases.

Figure 9 presents the flow rate estimations using the WarpCom method and the experimental setup and procedure that have been presented previously. If we compare the results obtained by estimation with the reference ones, we notice that the estimation performance stays almost unchanged for a flow rate between 1 and 3.4 m 3 /h. Above this value, there are more and more errors, having as main cause the vibrations of the pipe, that are present for a flow rate that is greater than 3 m 3 /h. We believe that these errors are not due to the methodology, but rather to the experimental physical setup, that allows that these vibrations cause vibrations on the transducers, leading to the observed errors in estimating the propagation times.

V. CONCLUSION

The results obtained by estimating the propagation times using the WarpCom method, correlated with the physics of the experimental setup, allow good estimations for the hydraulic parameters, as the flow rate and the temperature of the water. The results we obtained for real-world signals show the feasibility of the principle of estimating the hydraulic parameters by measuring the differences of propagation times.

Figure 1 .

 1 Figure 1. Illustration of the propagation of an acoustic wave

Figure 2 .

 2 Figure 2. The dependence of the speed of sound with temperature

Figure 3 .

 3 Figure 3. The flow rate measurement principle

Figure 4 .

 4 Figure 4. Experimental setup for measuring the temperature

Figure 5 .

 5 Figure 5. Experimental setup for measuring the flow rate

Figure 6 .

 6 Figure 6. Block diagram of the WarpCom procedure

Figure 7 .

 7 Figure 7. Polynomial modeling for a pulse train

  ) and estimated (red) sound speed

Figure 8 .

 8 Figure 8. Results obtained for the temperature

Figure 9 .

 9 Figure 9. Results obtained for the flow rate