Yacine Kessaci
email: yacine.kessaci@inria.fr

Mohand Mezmaz
email: mohand.mezmaz@umons.ac.be

Nouredine Melab
email: nouredine.melab@lifl.fr

E-G Talbi
email: talbi@lifl.fr

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

El-Ghazali Talbi

Parallel Evolutionary Algorithms for Energy Aware Scheduling

Reducing energy consumption is an increasingly important issue in computing and embedded systems. In computing systems, minimizing energy consumption can significantly reduces the amount of energy bills. The demand for computing systems steadily increases and the cost of energy continues to rise. In embedded systems, reducing the use of energy allows to extend the autonomy of these systems. In addition, the reduction of energy decreases greenhouse gas emissions. Therefore, many researches are carried out to develop new methods in order to consume less energy. This chapter gives an overview of the main methods used to reduce the energy consumption in computing and embedded systems.

As a use case and to give an example of a method, the chapter describes our new parallel bi-objective hybrid genetic algorithm that takes into account the completion time and the 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.

Introduction

Computers use a significant and growing portion of the energy in the world. Therefore, energy-aware computing is crucial for large-scale systems that consume considerable amount of energy and embedded systems that utilize battery for their power. A recent study on power consumption by servers [START_REF] Koomey | Estimating total power consumption by servers in the u.s. 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].

In this chapter, we present some important works done in the literature to reduce energy consumption. These works are classified according to three criteria. The first criterion concerns the optimization method used to minimize the consumed energy. This method can be mono-objective or multi-objective approach. The second criterion concerns the level of the system on which an approach is based. Indeed, a method can be based on the hardware or the software part of a system. The third and last criterion concerns the type of system to which the approach is intended to be used. An approach can be developed for computing systems or embedded systems.

In order to give an example of a method and as a use case, we describe our 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] Young | Minimizing energy consumption for precedenceconstrained 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] Topcuouglu | 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 combination 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.

The remainder of the chapter is organized as follows. Section 2 gives an overview of our classification. Section 3 explains the classification according to the optimization method used to minimize the consumed energy. Section 4 describes the classification according to the level of the system on which an approach is based. Section 5 explains the classification according to the type of system to which the approach is intended to be used. Section 6 presents the application, system, energy and schedul-ing models used in our use case. Our algorithm is presented in Section 7. The results of our comparative experimental study are discussed in Section 8. The conclusion is drawn in Section 9.

Energy aware approaches

This section presents an overview of the main methods published in the literature to reduce energy consumption. Our state of the art does not address the techniques used by computer system manufacturers. Indeed, manufacturers use some techniques in the design of electronic circuits and components to reduce their energy consumption. Our classification only focuses on methods intended to be exploited by the users of computing and embedded systems. As shown in Fig. 1, these methods can be classified according to three criteria.

Fig. 1 The classification criteria used

The first criterion concerns the type of the used optimization. This optimization can be mono-objective or multi-objective. In multi-objective optimization, the main approaches found in the literature are based on an aggregation, a lexicographic or a Pareto fitness function.

The system level is the second criterion on which methods can be classified. To reduce energy consumption, the methods are based either on software or hardware part of a system. In the hardware level, the techniques used are generally virualisation and consolidation based approaches. While in the software level, two techniques can be identified to address the task scheduling problem: Static methods where scheduling is done before the execution of a program, and dynamic methods, where the appropriate scheduling is calculated during the execution of a program.

The third and last criterion used in our classification is the type of system to which the approach is intended to be used. In this criterion, a system can be either a computing or an embedded system. Reducing energy consumption in embedded systems aims to increase the autonomy of devices. While the goal of reducing energy consumption in computing systems is to decrease the cost of the energy and the greenhouse gas emissions. Table 2 gives some examples of approaches with their classification using our taxonomy.

Method

Level System Optimization DPM&DVS [START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF] Hardware(Dyn.) E.S. Mono-objective DynamicAssgn [START_REF] Kang | Energy-efficient dynamic scheduling on parallel machines[END_REF] Hardware(Dyn.) C.S.(distributed comp.) Multi-objective(Lex.) ECS [START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF] Hardware(Dyn.) C.S.(high-performance comp.) Multi-objective(Agg.) ECTC and MaxUtil [START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF] Software(Cons.) C.S.(cloud comp.) Mono-objective Jitter [START_REF] Vincent | Just-intime dynamic voltage scaling: Exploiting inter-node slack to save energy in mpi programs[END_REF] Hardware(Dyn.) C.S.(high-performance comp.) qyut Multi-objective(Lex.) MMF-DVFS [START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF] Hardware(Stat.) C.S.(distributed comp.) Multi-objective(Lex.) MO [START_REF] Berl | A virtualized energy-efficient office environment[END_REF] Software(Virt.) C.S.(cloud comp.) Multi-objective(Lex.) PSAGA/PGA [START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF] Hardware(Dyn.) E.S. Mono-objective TDVAS [START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF] Hardware(Stat.) C.S.(cluster comp.) Multi-objective(Lex.)

Table 1 Examples of approaches with their classification

Optimization approaches

An approach used to optimize the energy consumption can be mono-objective or multiobjective. The mono-objective approaches, such as [START_REF] Hlavacs | Distributed energy efficiency in future home environments[END_REF][START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF][START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF][START_REF] Liu | Distributed energy-efficient scheduling for data-intensive applications with deadline constraints on data grids[END_REF][START_REF] Nathuji | Virtualpower: Coordinated power management in virtualized enterprise systems[END_REF][START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF][START_REF] Zhuo | Energy-efficient dynamic task scheduling algorithms for dvs systems[END_REF], only aim to optimize the energy consumption with assuming some constraints, especially on the completion time. In addition to the consumed energy, multi-objective methods are also trying to optimize other objectives such as quality of service (QoS).

QoS is usually expressed in terms of makespan. The energy consumption depends heavily on the QoS offered. Indeed, the minimization of energy consumption and the maximizing the QoS are conflicting objectives. The goal of a multi-objective optimization method is to find a good compromise between these two objectives. In addition, the cost of a service can be another objective to consider. More the QoS and energy consumption increase, more their cost increases. Table 3 gives some examples of the optimized objectives in different methods found in literature.

Method Optimization Criteria DPM&DVS [START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF] Mono-objective Energy DynamicAssgn [START_REF] Kang | Energy-efficient dynamic scheduling on parallel machines[END_REF] Multi-objective(Lex.) Energy and completion time ECS [START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF] Multi-objective(Agg.) Energy and completion time ECTC and MaxUtil [START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF] Mono-objective Energy Jitter [START_REF] Vincent | Just-intime dynamic voltage scaling: Exploiting inter-node slack to save energy in mpi programs[END_REF] Multi-objective(Lex.) Energy and completion time MMF-DVFS [START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF] Multi-objective(Lex.) Energy and completion time MO [START_REF] Berl | A virtualized energy-efficient office environment[END_REF] Multi-objective(Lex.) Energy and QoS PSAGA/PGA [START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF] Mono-objective Energy TDVAS [START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF] Multi-objective(Lex.) Energy and completion time As shown in Fig. 2, it is possible to identify three main categories according to the multi-objective approach chosen: Aggregation, lexicographic and Pareto approaches.

Aggregation approach

The aggregation (or weighted) method is one of the first and most used methods for generation of Pareto optimal solutions. It consists in using an aggregation function to transform a multi-objective problem into a mono-objective problem by combining the various objective functions into a single objective function generally in a linear way. The obtained results in the resolution of the problem depend strongly on the parameters chosen for the weight vector.

In [START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF], the authors address the task scheduling problem on heterogeneous computing systems (HCSs) and propose an energy-conscious scheduling heuristic (ECS) that takes into account the completion time and energy consumption. The heuristic proposed in [START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF] tries to balance these two performance goals using a novel objective function called relative superiority (RS).

Lexicographic approach

In this traditional approach, the search is carried out according to a given preference order of the objectives. This order defines the significance level of the objectives. Then, a set of mono-objective problems are solved in a sequential manner. If the problem associated with the most significant objective function has a unique solution, the search provides the optimal solution and stops. Otherwise, the problem associated with the second most significant objective function is solved. The same stopping criteria and process are iterated until the treatment of the last function.

In [START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF], the authors investigate the task scheduling problem, and propose a heuristic called maximum minimum frequency DVFS (MMF-DVFS). The proposed method operates in two stages. The goal of the first stage is to find a schedule of tasks that minimizes the makespan. The second stage tries to find the right setting of the processor to minimize energy consumption without changing the makespan of the schedule provided by the first stage.

Pareto approach

The Pareto approaches use the concept of dominance in the fitness assignment, contrary to the other approaches that use a scalarization function or treat the various objectives separately. The main advantage of the Pareto approaches is that they do not need the transformation of the multi-objective problem into a mono-objective problem. In a single run, they are able to generate a diverse set of Pareto solutions in the concave portions of the convex hull of feasible objective space.

For example, [START_REF] Miao | Energy saving task scheduling for heterogeneous cmp system based on multi-objective fuzzy genetic algorithm[END_REF] and [START_REF] Srikantaiah | Energy Aware Consolidation for Cloud Computing[END_REF] are one of the methods that use a Pareto approach. In [START_REF] Miao | Energy saving task scheduling for heterogeneous cmp system based on multi-objective fuzzy genetic algorithm[END_REF], the authors propose a multi-objective fuzzy genetic algorithm to optimize the energy saving scheduling tasks on heterogeneous chip multi-processor (CMP) system. In this algorithm, the Pareto set includes chromosomes with the shortest task execution time and the lowest system energy. In [START_REF] Srikantaiah | Energy Aware Consolidation for Cloud Computing[END_REF], the authors study the inter-relationships between energy consumption, resource utilization, and processor utilization.

System levels

To reduce energy consumption, various issues such as resource management in both software and hardware must be addressed. The system level (i.e. hardware or software) used by a method is another criterion on which methods can be classified. Software approaches are mainly based on virtualization or task consolidation. Virtualization consists to run on a single computer multiple operating systems as if they are running on separate computers. Hardware approaches use the opportunity offered by manufacturers in modern processors to adjust the voltage and frequency. This adjustment can vary the processor performance and thus its energy cost.

Hardware level

Dynamic voltage frequency scaling (DVS) is a powerful management technique in computer architecture. This technique enables processors to dynamically increase or decrease voltage supply levels (VSLs). DVS to decrease voltage is known as under-volting, and DVS to increase voltage is known as over-volting. Under-volting is done in order to conserve power, and over-volting is done in order to increase computer performance. Dynamic frequency scaling (DFS) is another powerful conservation technique that works on the same principles as DVS. This technique reduces the number of instructions a processor can issue in a given amount of time, thus reducing performance. Voltage and frequency scaling are often used together to save power. When used in this way it is commonly known as dynamic voltage and frequency scaling (DVFS). Many methods, such as [START_REF] De | Trade-offs between voltage scaling and processor shutdown for low-energy embedded multiprocessors[END_REF][START_REF] Vincent | Just-intime dynamic voltage scaling: Exploiting inter-node slack to save energy in mpi programs[END_REF][START_REF] Kang | Energy-efficient dynamic scheduling on parallel machines[END_REF][START_REF] Kimura | Emprical study on reducing energy of parallel programs using slack reclamation by dvfs in a power-scalable high performance cluster[END_REF][START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF][START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF][START_REF] Miao | Energy saving task scheduling for heterogeneous cmp system based on multi-objective fuzzy genetic algorithm[END_REF][START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF][START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF][START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF][START_REF] Zhuo | Energy-efficient dynamic task scheduling algorithms for dvs systems[END_REF], try to exploit this technique to reduce energy consumption.

In [START_REF] Vincent | Just-intime dynamic voltage scaling: Exploiting inter-node slack to save energy in mpi programs[END_REF], the authors present a system called Jitter, which reduces the frequency on nodes that are assigned less computation and therefore have slack time. This saves energy on these nodes, and the goal of Jitter is to attempt to ensure that they arrive just in time so that they avoid increasing overall execution time. [START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF] proposes an energy-efficient scheduling algorithm (TDVAS) using the dynamic voltage scaling technique to provide significant energy savings for clusters. The TDVAS algorithm aims at judiciously leveraging processor idle times to lower processor voltages (i.e., the dynamic voltage scaling technique or DVS), thereby reducing energy consumption experienced by parallel applications running on clusters. Reducing processor voltages, however, can inevitably lead to increased execution times of parallel task. The salient feature of the TDVAS algorithm is to tackle this problem by exploiting tasks precedence constraints.

In [START_REF] Kimura | Emprical study on reducing energy of parallel programs using slack reclamation by dvfs in a power-scalable high performance cluster[END_REF], the authors proposed a new algorithm that reduces energy consumption in a parallel program executed on a power-scalable cluster using DVFS. Whenever the computational load is not balanced, parallel programs encounter slack time, that is, they must wait for synchronization of the tasks. This algorithm reclaims slack time by changing the voltage and frequency, which allows a reduction in energy consumption without impacting on the performance of the program.

[10] also presents a novel dynamic algorithm for remapping tasks for energy efficient scheduling of DAG based applications for DVS enabled systems.

There are a lot of work on task scheduling problem using the dynamic voltage and frequency scaling technique in heterogeneous computing and real-time embedded systems. As previously mentioned, the main idea of the methods proposed in the literature is to change the voltage level to reduce energy consumption. These voltage level changes are made by taking into account timing information of processors, such as idle or slack time, and timing information of tasks, such as task deadline, task release time and task execution time. According to the availability of timing information of tasks, it is possible to distinguish two categories of methods to solve an energy-aware task scheduling problem: Static scheduling and dynamic scheduling.

In static scheduling [START_REF] De | Trade-offs between voltage scaling and processor shutdown for low-energy embedded multiprocessors[END_REF][START_REF] Miao | Energy saving task scheduling for heterogeneous cmp system based on multi-objective fuzzy genetic algorithm[END_REF][START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF][START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF][START_REF] Zhuo | Energy-efficient dynamic task scheduling algorithms for dvs systems[END_REF], timing information of tasks and timing information of processors are available in compile time before the deployment of tasks. This static scheduling involves most of the large-scale computational problems such as object recognition in machine vision applications, chemistry and bioinformatics. Having this information allows static schedulers to be developed by maximizing processor utilization with meeting all timing information of tasks.

In dynamic scheduling [START_REF] Vincent | Just-intime dynamic voltage scaling: Exploiting inter-node slack to save energy in mpi programs[END_REF][START_REF] Kang | Energy-efficient dynamic scheduling on parallel machines[END_REF][START_REF] Kimura | Emprical study on reducing energy of parallel programs using slack reclamation by dvfs in a power-scalable high performance cluster[END_REF][START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF][START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF][START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF], also called online or real-time scheduling, the deadlines of tasks are available in compile time but their release and execution times should be estimated during run time. This dynamic scheduling not only involves dynamic large-scale approximation and optimization such as in weather forecasting and search algorithms but also is used in most of power-aware devices like laptops, wireless sensors and cell phones.

Software level

In software level, two main techniques are used to reduce energy consumption: Virtualization and task consolidation. Operating system-level virtualization is a server method where the operating system allows for multiple isolated user-space instances. Virtualization can be seen as splitting an underlying hardware entity into smaller identical virtual entities running isolated from each other. Virtualization allows to reduce energy consumption

In [START_REF] Hlavacs | Distributed energy efficiency in future home environments[END_REF], an architecture for virtualizing and sharing hardware resources in future home environments is presented. The architecture aims at utilizing existing home resources in such a way that the consumed energy is minimized and the energy is efficiently used. A fully decentralized management system is proposed, interconnecting possibly thousands of homes in a peer-to-peer like manner. Energy optimization is done in a decentralized way by converging to a global energy optimum based on energy and performance metrics that have been defined. [START_REF] Berl | A virtualized energy-efficient office environment[END_REF] investigates the energy consumption in office environments and discusses the potential of energy savings. An energy-efficient office management approach is suggested, based on resource virtualization, power management, and resource sharing.

Like virtualization, task consolidation is an effective method to increase resource utilization and in turn reduces energy consumption. The task consolidation problem is the process of assigning a set N of n tasks (service requests or simply services) to a set R of r resources without violating time constraints. This technique aims to maximize resource utilization by minimizing energy consumption. Indeed, a resource allocation strategy that takes into account resource utilization would lead to better energy efficiency.

[24] exposes some of the complexities in performing consolidation for power optimization, and proposes some research directions to address the challenges involved. In [START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF], the authors present two energy-conscious task consolidation heuristics. These two heuristics aim to maximize resource utilization and explicitly take into account both active and idle energy consumption. The proposed heuristics as-sign each task to the resource on which the energy consumption for executing the task is explicitly or implicitly minimized without the performance degradation of that task.

System types

As shown in Fig. 4, two types of systems can be identified: Embedded systems and computing systems. The issue of reducing energy consumption interests these two types of systems for different reasons. In embedded systems, the main goal of reducing energy consumption is to increase the autonomy of devices, while in computing systems, the main goal is to reduce the energy cost.

Embedded systems

An embedded system can be defined as an electronic and autonomous system which is dedicated to a specific task. The resources of an embedded system are limited. This limitation is usually spatial (limited size) and energy (consumption restricted). Many of today's embedded systems, such as wireless and portable devices rely heavily on the limited power supply. Therefore, energy efficiency becomes one of the major design concerns for embedded systems. Many methods, such as [START_REF] De | Trade-offs between voltage scaling and processor shutdown for low-energy embedded multiprocessors[END_REF][START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF][START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF][START_REF] Zhuo | Energy-efficient dynamic task scheduling algorithms for dvs systems[END_REF], are developed to reduce the energy consumption of embedded systems.

In [START_REF] De | Trade-offs between voltage scaling and processor shutdown for low-energy embedded multiprocessors[END_REF], some scheduling heuristics are presented that determine the best trade-off between three techniques: DVS, processor shutdown, and finding the optimal number of processors. Indeed, when peak performance is unnecessary, DVS can be used to reduce the dynamic power consumption of embedded multiprocessors. However, static power consumption is expected to increase significantly. Then it will be more effective to limit the number of employed processors, and use a combination of DVS and processor shutdown. [START_REF] Simunic | Dynamic voltage scaling and power management for portable systems[END_REF] presents an optimal approach based on the time-indexed semi-Maxkov decision processes (TISMDP) for optimizing power management policies of embedded systems.

In [START_REF] Lin | Parallel genetic algorithms for dvs scheduling of distributed embedded systems[END_REF], the authors developed a genetic algorithm for the DVS scheduling problem. They also describe a Parallel genetic algorithm for finding better schedules with less time by parallelizing the genetic algorithm. A hybrid parallel algorithm is also developed to further improve the search ability of parallel genetic algorithm by combining the parallel genetic algorithm with the technique of simulated annealing.

Computing systems

Demand for cloud computing, high performance computing, cluster computing, etc. is growing. The industry is responding to the demand by more powerful systems, and consequently, creating highly energy-intensive products. The energy consumption of these systems is associated with various environmental, system performance and monetary issues.

Scheduling parallel applications on computing systems, especially on large-scale computing systems, is challenging due to high energy consumption and significant communication latencies. Therefore, conserving energy consumption and reducing schedule lengths are two important concerns in the design of economical and environmentally friendly systems. It is only recently that much attention has been paid to energy consumption in scheduling. [START_REF] Kang | Energy-efficient dynamic scheduling on parallel machines[END_REF] and [START_REF] Babaii Rizvandi | Linear combinations of dvfs-enabled processor frequencies to modify the energy-aware scheduling algorithms. Cluster Computing and the Grid[END_REF] addresse the task scheduling problem on distributed computing systems. [START_REF] Young | Minimizing energy consumption for precedenceconstrained applications using dynamic voltage scaling[END_REF] and [START_REF] Miao | Energy saving task scheduling for heterogeneous cmp system based on multi-objective fuzzy genetic algorithm[END_REF] investigate the same problem on respectively heterogeneous computing systems and chip multi-processor systems. With the chip multi-processor being more and more widespread used in the laptop, desktop and data center area, the power-performance scheduling issues are becoming challenges. In [START_REF] Ruan | An energy-efficient scheduling algorithm using dynamic voltage scaling for parallel applications on clusters[END_REF], the authors address the scheduling parallel applications on large scale clusters. In the past decade cluster computing platforms have been widely applied to support a variety of scientific applications. [START_REF] Liu | Distributed energy-efficient scheduling for data-intensive applications with deadline constraints on data grids[END_REF] proposes an energy efficient algorithm to schedule real-time tasks with data access requirements on grids. Taking into account both data locations and application properties, the authors design a distributed energy-efficient scheduler that aims to seamlessly integrate the process of scheduling tasks with data placement strategies to provide energy savings.

Consolidation of applications in cloud computing environments presents an important approach to streamline resource usage and improve energy efficiency. Based on the fact that resource utilization directly relates to energy consumption, the authors in [START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF] have modeled their relationship and developed two energy-conscious task consolidation heuristics. [START_REF] Srikantaiah | Energy Aware Consolidation for Cloud Computing[END_REF] outlines the challenges in finding effective solutions to the consolidation problem on cloud computing systems. In [START_REF] Hlavacs | Distributed energy efficiency in future home environments[END_REF], an architecture for sharing computing resources among home environments in a peer-to-peer manner is proposed in order to improve the energy efficiency. [START_REF] Nathuji | Virtualpower: Coordinated power management in virtualized enterprise systems[END_REF] explores how to integrate power management mechanisms and policies with the virtualization tech-nologies being actively deployed in large-scale datacenters to address costs and limitations in cooling or power delivery.

Problem modeling

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

System model

The target system used in this work consists of a set P of p heterogeneous processors/machines that are fully interconnected. 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 3). Since clock frequency transition overheads take a negligible amount of time (e.g., 10µs-150µs [9], [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 i f n i = n entry EFT (MIP(n i ), p k ) + c MIP(n i ),i otherwise (1) EFT (n i , p j ) = EST (n i , p j ) + w i, j (2)
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 AFT (n i , p j ) can be different from its earliest start and finish times, EST (n i , p j ) and EFT (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 Fig. 5 with its details in Table 4 and Table 5. The values presented in Table 5 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, shortcircuit 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 ( 3
)
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 (3) 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 ( 4
)
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{AFT (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.

A parallel evolutionary algorithm

In this section, our new parallel bi-objective hybrid approach is presented.

Hybrid approach

In our approach illustrated in Fig. 6, 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 ).

Fig. 6 Our hybrid GA (GA and ECS)

The new approach we propose is based on ECS which is not a population-based heuristic. ECS tries to construct in a greedy way one solution using three components.

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

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

• And a 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 6. 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 Fig. 7). Fig. 7 The mutation operator As illustrated in Fig. 8, 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.

Fig. 8 The crossover operator

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 Fig. 9, 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 Fig. 10, our parallelization is based on the deployment of the approach using Fig. 9 The cooperative island approach 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. The new approach is experimented on about 10,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.).

Experiments have been performed on a grid of three clusters. 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 6. Each instance is solved twice. The first resolution is done with ECS, and the 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.

Experiments show that our approach improves on average the results obtained by ECS. Indeed, the energy consumption is reduced by 47.49% and the makespan reduced by of 12.05%. In addition, our experiments show 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 experiments, presented in this section, focus only on the large instances of Table 6. 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. Fig. 11 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 Fig. 11, 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 Fig. 11, 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 8, 32, or 64, and the heterogeneity of processors is 100, 200 or random. The population of the GA contains 20 chromosomes. Therefore, 21 computing cores are used to solve each instance (20 cores to run the ECSs and 1 core to run the GA). In our case, an experiment can not have a speedup greater than 21.

The average speedup obtained is 13.06. Our experiments show that the speedup increases proportionally to the number of processors on which the precedenceconstrained parallel application is run, and the CCR and the heterogeneity of processors do not impact significantly the quality of the acceleration of our approach.

Conclusions

In this chapter, we presented a classification of different methods used in literature to reduce energy consumption. Our classification is made according to three criteria. These criteria are the optimization method used to minimize the consumed energy (i.e. mono-objective or multi-objective optimization), the level of the system on which an approach is based (i.e. hardware or software level), and the type of system to which the approach is intended to be used (i.e. computing or embedded systems).

As a use case and to give an example of a method, we presented a new parallel biobjective hybrid genetic algorithm to minimize energy consumption and makespan. The energy saving of our approach exploits the dynamic voltage scaling (DVS). According to our classification, the new method can be considered as an optimizing multi-objective method with a Pareto approach. It uses the hardware part of the system in a dynamic way. Our method is intended to be used in computing systems.

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] Topcuouglu | 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 npcompleteness[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.

Therefore, one of the main perspectives of the work presented in this chapter 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.

Fig. 2

 2 Fig. 2 Classification according to the optimization approaches

Fig. 3

 3 Fig. 3 Classification according to the system levels

Fig. 4

 4 Fig. 4 Classification according to the system types

Fig. 5 A

 5 Fig. 5 A simple task graph

Fig. 10

 10 Fig. 10 Illustration of the multi-start approach

Fig. 11 S

 11 Fig. 11 S-metric value according to the number of islands

Table 2

 2 Examples of the optimized objectives in different methods

Table 3

 3 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 4

 4 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
	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 5

 5 Task Priorities

Table 6

 6 Table 6 summarizes the parameters used in our experiments. 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

Acknowledgments

Experiments presented in this chapter 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.


Mitsuhiro Hayashibe
email: hayashibe@lirmm.fr

Gentiane Venture
email: venture@cc.tuat.ac.jp

Ko Ayusawa
email: ayusawa@ynl.t.u-tokyo.ac.jp

Yoshihiko Nakamura
email: nakamura@ynl.t.u-tokyo.ac.jp

Muscle Strength and Mass Distribution Identification Toward Subject-Specific Musculoskeletal Modeling

In current biomechanics approach, the assumptions are commonly used in body-segment parameters and muscle strength parameters due to the difficulty in accessing those subject-specific values. Especially in the rehabilitation and sports science where each subject can easily have quite different anthropometry and muscle condition due to disease, age or training history, it would be important to identify those parameters to take benefits correctly from the recent advances in computational musculoskeletal modeling. In this paper, Mass Distribution Identification to improve the joint torque estimation and Muscle Strength Identification to improve the muscle force estimation were performed combined with previously proposed methods in muscle tension optimization. This first result highlights that the reliable muscle force estimation could be extracted after these identifications. The proposed framework toward subject-specific musculoskeletal modeling would contribute to a patient-oriented computational rehabilitation.

I. INTRODUCTION

A. Background

An understanding of the whole body human motion is a complex process that requires information at different scales. The observation of a motion reflects in fact the brain activity, the neural system activity, and the musculo-skeletal system activity. To fully understand the mechanisms and the internal network of human motion generation, it is thus mandatory to understand each element and the connections between these elements in a multidisciplinary framework. Muscles are the basic elements to actuate the joints in the human body [START_REF] Mcmahon | Muscles, Reflexes, and Locomotion[END_REF]. Analyzing how they activate and generate forces for a specific motion would significantly improve several research areas ranging from physical therapy to neuro-rehabilitation. In addition, quantitative analysis between neural activity and muscle force can contribute to the design of Functional Electrical Stimulation (FES) for paralyzed muscles and myoelectrical limb control, and also to the design of robotic prosthetic limbs and exoskeletons.

The modeling of the human body and the computation of human dynamics is an active field in robotics and biomechanics. The general musculoskeletal model of whole body and its dynamics computation method were established and are now available as OpenSim software [START_REF] Delp | A computational framework for simulating and analyzing human and animal movement[END_REF] [START_REF] Delp | OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement[END_REF]. The efficient computational algorithms for musculoskeletal model were pursued for inverse and forward dynamics on the basis of efficient multibody dynamics computations [START_REF] Nakamura | Somatosensory Computation for Man-Machine Interface From Motion-Capture Data and Musculoskeletal Human Model[END_REF] [START_REF] Yamane | Estimation of physically and physiologically valid somatosensory information[END_REF].

In this paper, we discuss the muscle strength and mass distribution identification toward subject-specific musculoskeletal modeling. As mentioned above, many researches on musculoskeletal modeling have been reported especially in the biomechanics society. These techniques are getting mature to be applied for the clinical rehabilitation study and other specific purposes. However, to our knowledge, most of the works use the literature-based anthropometric index for segment mass distribution. The work of de Leva [START_REF] De Leva | Adjustments to Zatsiorsky-Seluyanov's segment inertia parameters[END_REF] which adjusted the Zatsiorsky model is commonly used for segment definitions to obtain "typical" mass distribution. The subject weight, height are used as the reference to compute the corresponding mass of each segment derived from anthropometric table. As it depends on the total mass, it should not be able to represent the subject-specific difference of mass distribution in the case where the subjects have the same weights and heights, the one is an athlete, and the other is an elderly person. In addition, regarding muscle isometric force table, the work of Delp [START_REF] Delp | OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement[END_REF] is quite famous and often used when the muscular force is computed based on electromyography (EMG). This musculotendon parameters database is meaningful as it is normally difficult to access these typical values. However, we can not assume that every subject would have similar muscle strength especially in the rehabilitation and sports science where each subject has quite different muscle condition due to disease, age or training history. In fact, the group even in the same generation and healthy people can easily have different muscle strength depending on daily exercises and the preference of sports.

B. Related works

Inverse dynamics of musculoskeletal models is commonly used to estimate the joint torques and muscle tensions from motion, ground contact force, and/or EMG measurements. Known problems with the inverse dynamics approach is that one must know the inertia and mass of each body segment in order to estimate joint moments correctly [START_REF] Buchanan | Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command[END_REF]. Typically these are estimated using values from cadavers and scaled using simplistic scaling rules, the accuracies of which are rarely verified. This scaling technique is employed also in OpenSim [START_REF] Delp | OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement[END_REF]. A novel identification method was proposed and implemented for the human body-segment mass-parameters [START_REF] Venture | Real-time identification and visualization of human segment parameters[END_REF]. This method can contribute to reduce errors in the inverse dynamics result significantly. Its realtime implementation is also available [START_REF] Ayusawa | Real-time implementation of physically consistent identification of human body segments[END_REF].

The resultant joint torques from inverse dynamics are net values contributed from every muscle. As the numbers of muscles span for each joint is greater than the degrees of freedom for each joint, the estimation of each muscle force is a redundant problem. Thus, usage of EMG has been studied by many researchers in biomechanics. EMG-based model relies on measured muscle activity to estimate muscle force. Potvin et al. [START_REF] Potvin | gMechanically corrected EMG for the continuous estimation of erector spinae muscle loading during repetivtive lifting[END_REF] derived the relationship between the EMG signal and muscle force to estimate the muscle load during lifting motion. EMG-based models have been used to estimate torques around the knee joint by calibrating the muscle parameters for individual subjects [START_REF] Lloyd | An emg-driven musculoskeletal model to estimate muscle forces and knee joint moment in vivo[END_REF]. However, previously proposed methods for the estimation of forces based on EMG with calibration were about only one DOF as it requires dynamometer to obtain reliable joint torque. Recently, EMG-driven musculoskeletal model was proposed for multiple DOFs of the joints [START_REF] Sartori | An EMG-driven Musculoskeletal Model of the Human Lower Limb for the Estimation of Muscle Forces and Moments at the Hip, Knee and Ankle Joints in vivo[END_REF]. This method uses Opensim to obtain joint torques from motion capture data.

Drawback of complete EMG-driven musculoskeletal model is that it is unrealistic to access all the inner muscles. In addition, the tensions computed from EMG data are not reliable due to the noise in EMG measurements and uncertain muscle model parameters. Thus, Yamane and Nakamura proposed to obtain physically and physiologically reasonable whole-body muscle tensions by forming an optimization problem such that the error of the tension-torque transformation equation becomes minimum while respecting the physiological muscle dynamics and the muscle activity when it is available from EMG [START_REF] Yamane | Estimation of physically and physiologically valid somatosensory information[END_REF]. All inner muscles are included in the optimization. However, subject-specific muscle strength is not considered yet. Identification method to identify subject-specific muscle states and maximal force parameter was developed for FES in [START_REF] Hayashibe | Nonlinear Identification Method Corresponding to Muscle Property Variation in FES -Experiments in Paraplegic Patients[END_REF]. In biomechanics, normally dynamometer is used for subject-specific muscle strength information.

This work is the first trial to make the identification on subject-specific muscle strength and mass distribution based on motion capture. The concept of this paper is summarized in Fig. 1. The white colored parts are the standard processes in the current musculoskeletal analysis. As the subjectspecific mass distribution is identified using kinematic information and ground reaction force, the computed joint torques from the inverse dynamics of musculoskeletal models can be improved to output reliable physically consistent torques. In addition, as the subject-specific muscle strength is identified with the above torque and EMG signal, the optimization result of muscle tensions can be improved to output reliable physiologically consistent torques reflecting subject-specificity.

II. MASS-DISTRIBUTION IDENTIFICATION

The equations of motions of a moving system can be written as Eq. 1 [START_REF] Yoshida | Moving base robotics and reaction management control[END_REF] kinematic chains constituting the whole system. The baselink can be chosen arbitrarily, however for practical reasons it is convenient to select the lower torso. The chains are formed with the limbs, the upper-torso, the head.

[

H oo H oc H co H cc ] [ qo θ ] + [ b o b c ] = [ 0 τ ] + Nc ∑ k=1 [ K ko K kc ] F k (1)
where:

• H ij (i, j = o, c) is the inertia matrix,
• q o is the vector of generalized coordinates which represents the 6 DOF of the base-link, • θ is the vector joint angles, • b i is the bias force vector including centrifugal, Coriolis and gravity forces, • τ is the vector of joint torques,

• N c is the number of contact point with the environment, • F k is the k th vector of external forces, • K ko and K kc are matrices which map F k to the generalized force vector.

In the aim of identification of the dynamics it is common to rewrite Eq.1 by separating the time-constant inertial parameters ϕ from the time varying functions Y [START_REF] Khalil | Modeling, identification and control of robots[END_REF], resulting in the identification model given by Eq.2. We have shown in our previous work [START_REF] Venture | Real-time identification and visualization of human segment parameters[END_REF] that it is possible to solve the identification problem using solely the upper equation of Eq.2, thus without requiring the joint torque information. This method is particularly adequate for the measurement of the human body dynamics.

[ Y o Y c ] ϕ = [ 0 τ ] + Nc ∑ k=1 [ K ko K kc ] F k (2)
where:

• Y = [ Y o Y c
] is the regressor, which is function of the system joint angles θ, velocity θ and acceleration θ, and of the vector of generalized coordinates q o and its derivatives; • ϕ is the vector of inertial parameters to estimate such that:

ϕ = [ ϕ T 0 ϕ T 1 •• ϕ T N ] T ( 3
)
• ϕ j is the vector of standard parameters for each link B j (j = 0 to N ), such that:

ϕ j = [ m i ms i,x ms i,y ms i,z I i,xx I i,yy I i,zz I i,yz I i,zx I i,xy ] T (4) • m j is the mass, • I j,xx , I j,
yy , I j,zz , I j,yz , I j,zx , I j,xy are the 6 independent components of the inertia matrix I j , • ms i,x , ms i,y , ms i,z are the first moments components of the vector ms j Because the measuring environment is not free from noise and measuring errors, and the modeling of the human body is an approximation, the identification results are biased. This bias, even with sufficiently exciting trajectories may lead to parameters that are not physically consistent. Namely, inertia matrices should be definite positive, masses should positives, and the center of a mass of a link appropriately located in the link. In order to guarantee the physical consistency of each parameters we have also included a constraint based on a simplified shape of the link [START_REF] Ayusawa | Real-time implementation of physically consistent identification of human body segments[END_REF]. For the geometric identification: lengths of limbs, center of joints, we use the marker positions captured from a T-pose to generate automatically the adequate model.

The model used for segment identification and the musculo-skeletal model differ in the number of degrees of freedom. This is because for identification we must consider the identifiability problem, and for the musculoskeletal we must take into account via-points of muscles. In the musculo-skeletal model several virtual links or fixed DOF are utilized, and the trunk is divided in links of small size. The identified parameters must then be adjusted to this model. The limbs, the head are globally the same, so we use directly the identified parameters. For the other links, we use a proportional distribution of mass and inertia to reflect the results of identification, while preserving the high number of DOF of the musculo-skeletal model.

III. TORQUE TO FORCE MAPPING IN IDENTIFIED MODEL

From the lower part of Eq. 1, once subject-specific inertia matrix is identified, more reliable joint torques τ can be computed with external forces using Newton-Euler inverse dynamics as in Eq.5. This process is indicated as ID in Fig. 1.

τ = [ H co H cc ] θ + b c - Nc ∑ k=1 K kc F k (5)
To obtain the relationship between the joint torques τ and the muscle element tensions f , the Jacobian matrix of the element lengths with respect to the generalized coordinates J ∈ R nmus×n dof is computed as Eq. 6.

J = ∂l ∂θ G ( 6
)
where:

• l ∈ R nmus is the muscle lengths, • θ G ∈ R n dof is the generalized coordinate,
• n mus and n dof are the number of muscles and DOF of the model. The method for computing J is detailed in [START_REF] Nakamura | Somatosensory Computation for Man-Machine Interface From Motion-Capture Data and Musculoskeletal Human Model[END_REF]. The principle of virtual works yield the following equation relating the muscle forces and the joint torques.

τ = J T f (7)
Solving Eq.7 for muscle forces f using joint torques τ is a highly redundant problem. Because the number of muscles n mus is greater than the total DOF n dof . The skeletal model used in this work has 83 DOF with 314 muscles.

IV. EMG PROCESSING AND MUSCLE CONTRACTION

DYNAMICS

The transformation from EMG to muscle activation is an essential process in Hill-type muscle models [START_REF] Hill | The heat of shortening and the dynamic constants in muscle[END_REF][12], it is dominant process because the estimated muscle force is assumed to be proportional to the muscle activation. One can refer to the detailed steps in [START_REF] Yamane | Estimation of physically and physiologically valid somatosensory information[END_REF]. First, mean of the raw EMG was taken to offset the baseline of EMG signal. The integrated EMG signal, called IEMG, normalized with the value of Maximum Voluntary Contraction (MVC) gives the activity level u. The muscle activity a is obtained by the following first-order differential equation:

ȧ = u -a T T = T a (u ≥ a) T = T d (u < a) ( 8
)
where T a and T d are the time constants for activation and deactivation.

The Hill-type model is used to estimate the reference muscle force f * with general form of the function by

f * = af l (l)f v ( l)F iso ( 9
)
where:

• F iso is the maximum isometric muscle force,

• l is the muscle length,

• f l (l) and f v ( l) gives the normalized force-length and force-velocity relationship respectively [START_REF] Buchanan | Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command[END_REF]. The force length relationship shows a Gaussian distribution around the optimal length l 0 and is formulated as

f l (l) = exp { - ( l -l 0 K l ) 2 } ( 10
)
where K l is a constant parameter. f v ( l) represents the relationship between velocity and normalized force. The muscle can contract at its maximum velocity v max without load and slows down as the load increases. This relationship is formulated as follows:

f v ( l) =        0 ( l ≤ -v max ) V sh (vmax+ l) V sh vmax- l (-v max ≤ l ≤ 0) V sh V shl vmax+V ml l) V sh V shl vmax+ l (0 ≤ l) (11)
where V sh ,V shl and V ml are constant parameters given by Stroeve [START_REF] Stroeve | Impedance Characteristics of a Neuro-Musculoskeletal Model of the Human Arm I: Posture Control[END_REF].

V. MUSCLE STRENGTH IDENTIFICATION

Considering the joint torque τ k which is generated by the muscle k,

τ k = J T k f k (12) = [-a k f l k (l k )f v k ( lk )J T k ]F isok (13) = C T k F isok ( 14
)
where J k ∈ R 1×n dof is the element Jacobian for muscle k. As f l k (l k ),f v k ( lk ),J T k can be obtained from the result of inverse kinematics and a k is obtained from EMG measurement for the muscles with EMG electrodes, C k is known for the muscles with EMG electrodes. We call C k as net contribution matrix.

As the human motion is generated by complex combination of multiple muscle forces, then it is not always the case but for certain tasks, particular joint torque among τ is expected being consisted of the contributions from particular muscles. In such case and with the assumption where synergist muscles exhibit similar activity levels like the muscle grouping in stretch reflex, the net contribution matrix C k are available if there is EMG measurement in the muscle group, certain elements in joint torque vector can be expressed as follows:

τ (j) = ∑ k=1 C T k (j)F isok (15) = ∑ Gemg C T k (j)F isok + d ( 16
)
where:

• (j) represents jth element of the vector,

• G emg is the EMG measured synergist muscle group,

• d is the constant residual term which takes into account the contribution from other elements. The subject-specific Muscle-Strength F isok is identified with linear regression.

After the identification of subject-specific Muscle-Strength F isok , a more reliable EMG-based reference muscle force f * can be computed. From the previous work [START_REF] Nakamura | Somatosensory Computation for Man-Machine Interface From Motion-Capture Data and Musculoskeletal Human Model[END_REF] about optimization criteria for the actuation redundancy problem in quadratic programming, we should now minimize the objective function to find a solution f keeping the error against EMG-based muscle tension which is physiologically consistent and the error against the resultant joint torque which is physically consistent and subject-specific. The optimization is formulated as follows:

Find f that minimizes

Z = 1 2 |τ -J T f | 2 + W emg 2 |f -f * | 2 (17)
subject to the following constraints:

E m f ≤ 0 (18)
where W emg is the weight for EMG-based force reliability coefficient and E m is a matrix which extracts the tensions of all elements from f . Defining the appropriate W emg is an issue involving the argument on how we can decide it, because normally f * is not a subject-specific estimation. For instance, when the literature based F iso is taken for the elderly subject or athlete, the right term of the objective function can lead the inconsistency in the sense of subject-specific muscle strength. Thus, the identification of muscle strength is a trial to make the EMG-based muscle tension subject-specific and physically consistent. Even with the same choice of the W emg , the solution f can be more reliable as its absolute value. In addition, the identified Muscle-Strength can be used for the performance index to evaluate the subject evolution in training and rehabilitation. In the optimistic scenario, software may automatically judge if the designated task can be performed by the patient who has particular motor function problem.

VI. EXPERIMENTAL RESULTS

A. Experimental setup

The human motions are recorded by a commercial optical motion capture system consisting of 10 cameras (Motion Analysis). 35 reflective optical markers pasted on the body are captured by these cameras. We employed modified Helen Hayes Hospital marker set. They are located at the defined anatomical rigid points to diminish the influence from nonrigid skin and muscle movement to insure accuracy of inverse kinematics computations. The contact forces are measured by two force-plates (Kistler). EMG data are recorded at 1kHz with a 16-channel wireless EMG system (DELSYS) that is synchronized with other measurements. Bipolar surface electrodes are placed on the following 16 muscles for both legs; Gluteus Maximus, Rectus Femoris, Vastus Lateralis, long head of Biceps Femoris, Semimembranosus, Tibialis Anterior, medial head of Gastrocnemius, Soleus as shown in Fig. 2.

First, subjects are requested to make Maximum Voluntary Contraction (MVC) for the extension and flexion in each joint of lower limbs. The corresponding EMG value was used as an initial reference for the normalization in EMG processing. Then, we record a T-pose of a few seconds for the geometric calibration. Next, a sequence of about 2 minutes of free exciting movements for the mass-distribution identification. For the better excitability for each joint DOF, subjects are requested to pay attention to articulate each body-segment. One session consisting of slow, normal and fast squat is repeated with different weighted loads (0kg, 4kg, 8kg) around the waist.

B. Identification of Mass-Distribution

The identified geometric model is used to compute the inverse kinematics. Each joint information and the global coordinates are then used to compute the regressor matrix and the transformation matrix for the forces. The contact forces are measured by two force plates directly. The model of the body used for identification is simplified compared to the musculo-skeletal model. It consists in 34 degrees of freedom located between 15 links. Except for the upper torso, the links correspond to actual segment, thus it is not a simplification for lower limbs. The upper-torso is considered as one link for a precise identification, rather than using a decomposition similar to the musculo-skeletal model in several links of extremely small size, and thus with little mass and inertia, therefore difficult to identify precisely. The global parameters of the upper torso are then redistributed on the small links to match the specificity of the subject. The identified results are visualized graphically as shown in Fig. 3. The green color indicates that the inertial parameters are well identified, based on statistic evaluation. The figure in the link is the identified mass in kilograms. The results are consistent with expected values: no aberrant value, no negative value, a good right/left balance. And there is a clear difference in the identified Mass-Distribution between candidate A and candidate B, which are of different morphologies.

C. Identification of Muscle-Strength

In this first trial of Muscle-Strength (M-S) identification, we focused on knee joint torque variation under squat motions. First, we needed to extract motions which have torque variation well correlated to EMG-measured synergist muscle groups to be applied in Eq.16. The example of muscle activation transition is shown in the up of Fig. 4. It was obtained from a squat motion with 8kg load. Correlation analysis in the level of the net contribution matrix was performed to judge which EMG-measured synergist muscle group have enough correlation and contribute to the torque production. The knee element of contribution matrix in Rectus Femoris, Vastus Lateralis had significant correlation Candidate B Candidate A Fig. 3. The results of the identification for the two candidates. The figure represents the identified mass of each link in kg unit, the green color shows link which parameters are identified with excellent accuracy, the blue color shows link which parameters are small. The red arrow shows the direct force plate measurement, the blue arrow is the reconstruction using the identified parameters. as shown in the bottom of Fig. 4. The squat motion with load made the situation where the knee extensor is required to be activated principally even in the knee bending phase. For Vastus Lateralis, the synergist muscle group: Vastus Medialis and Vastus Interioris are assumed to have same muscle activity. F iso was identified applying linear regression in Eq.16.

The estimation of muscle tensions was performed with the optimization as in Eq.17. The muscle force estimation with different weights W emg was carried out with original settings and M-S identified settings for F iso of Rectus Femoris and Vastus muscle group. The average error of left knee joint in the term of |τ -J T f | and the maximal estimated force in left Rectus Femoris of Candidate B during squat with 8kg were evaluated as in Table I. The maximal estimated force indicates the force scale solved in the optimization. F iso of Rectus Femoris with original model was 1kN . F iso was identified as 685N in Candidate B. F iso was identified as 802N in Candidate A. For reference, the parameter in Delp musculoskeletal table is 780N . Considering the difference of musculoskeletal geometrical path, it is well identified within similar range. Normally when W emg is increased, τ error can be increased because the solution needs to be found keeping physical consistency from ID and respecting EMG-based model output. τ error showed similar result both for the optimization with/without M-S identification. However, when M-S identification is not performed, the estimated force scale varied a lot depending on W emg . It is not preferable since the estimated force is influenced by the computation setting. It can be explained by the following: as the muscle reference force f * is not realistic, minimizing the right term in the objective function causes problem. Thus, the solution is easily influenced by the difference of W emg . In addition, as F iso is 1kN for original model, we can see the solution without M-S Identification is misled while increasing W emg . Thus, the estimated value is still usable as the comparative study with other muscles such as distribution ratio, but not for the individual analysis with different tasks. For very small W emg , even in M-S identified case, it started to give different solution, but it is normal as it starts not to respect to EMG-based model output. In any case, the optimal setting of W emg and W dyn is hard to be predicted before the analysis. Since the solution with M-S identification became stable with much larger range of W emg , this first results highlight the effectiveness of Muscle-Strength identification.

Next, the identified model was applied for the muscle forces estimation in another squat motions with 4kg load. Fig. 5 shows the estimated muscle forces of Rectus Femoris, Vastus Lateralis and Vastus Medialis with different speeds (first two series are normal speed, second two are slow and last two are fast). The corresponding visualizations of estimated muscle tensions at the indicated time instant are depicted in the bottom.

VII. CONCLUSION

In this paper we have presented a method that allows to reliable subject-specific muscle forces estimation with Muscle Strength and Mass Distribution Identification. The proposed method features:

• Mass Distribution Identification to improve the joint torque estimation accuracy, • Muscle Strength Identification to improve the muscle force estimation accuracy, • The above features result in subject-specific musculoskeletal dynamics computation to compensate the drawback of conventional computer aided analysis such as the assumptions in body-segment parameters and muscle strength parameters. The result as in Fig. 5 indicated the starting timing of contractions of knee extensor muscles is visualized in quite reasonable way and the estimated value is expected to be reliable since this model computation is considering the subject-specific mass distribution and subject-specific muscle strength of knee extensors along with EMG activities of surface muscles and all the inner muscles contributions in the optimization process.

However, in order to apply the proposed method to all the other joints, we still have open problem such as how we can automatically detect the motions which have high correlation between torque variation and net contribution matrix. The excitation index tracking algorithm in torquemuscle force mapping would be necessary in order to make this framework applicable in systematic way. The knowledge about modular control in synergist muscle group [START_REF] Ivanenko | Modular Control of Limb Movements during Human Locomotion[END_REF] is recently well accepted in neuroscience. It would help a lot to decrease the actuation redundancy problem. Muscle grouping to realize real-time computation is already implemented in [START_REF] Murai | Computationally Fast Estimation of Muscle Tension for Realtime Bio-feedback[END_REF]. Muscle grouping concept would help also for systematical detection of high correlated motion in torque-muscle force mapping.

Even though the above problems are required to be solved, the proposed framework toward subject-specific musculoskeletal modeling would contribute to a patientoriented computational rehabilitation. In addition to subjectspecificity, human characteristics are basically time variant, for instance, neuromuscular dynamics may vary according to muscle fatigue. In order to correspond to such time-varying characteristic, a model-based approach which is combined with adaptive identification would enable to bring a new modality in musculoskeletal modeling as in [START_REF] Zhang | FES-Induced Torque Prediction with Evoked EMG Sensing for Muscle Fatigue Tracking[END_REF].

Fig. 1 .

 1 Fig. 1. Scheme of Mass Distribution and Muscle Strength identification.

Fig. 2 .

 2 Fig. 2.Appearance of the experiment to measure EMGs and marker positions for lower limbs.

Fig. 4 .

 4 Fig.[START_REF] Nakamura | Somatosensory Computation for Man-Machine Interface From Motion-Capture Data and Musculoskeletal Human Model[END_REF]. The example of the muscle activation transition from a squat motion with 8kg load (up). The knee element transition of contribution matrix in Rectus Femoris, Vastus Lateralis during the same motion. They showed significant correlation with knee joint torque during the motion (bottom)

Fig. 5 .

 5 Fig. 5. Estimated muscle forces of Rectus Femoris, Vastus Lateralis and Vastus Medialis using the identified model in squat motion with 4kg load with different speeds (first two series are normal speed, second two are slow and last two are fast) (up). The corresponding visualization of estimated muscle tensions at the indicated time instant (bottom).

TABLE I KNEE

 I TORQUE ERROR AND ESTIMATED FORCE SCALE FOR RECTUS FEMORIS WITH DIFFERENT Wemg IN STATIC OPTIMIZATION error showed similar result both for the optimization with/without M-S identification. The estimated force is less influenced by Wemg with M-S identification.

		with M-S Ident	no M-S Ident
	Wemg	τ error max force	τ error max force
	1	1.1	632	0.87	940
	5 × 10 -1	0.99	626	0.78	938
	10 -1	0.59	595	0.49	871
	5 × 10 -2	0.45	591	0.39	832
	10 -2	0.26	585	0.23	543
	5 × 10 -3	0.23	580	0.20	433
	10 -3	0.19	483	0.17	337
	τ

VIII. ACKNOWLEDGMENTS

The authors acknowledge the use of sDIMS, multibody and musculoskeletal dynamics computation software [4] [5], developed in the University of Tokyo.

Authors acknowledge the support provided by the INRIA-JSPS joint program AYAME (2010-2012).

M.