Member, IEEE Lei Mo

Angeliki Kritikakou

Member, IEEE Olivier Sentieys

Transactions on Computer-Aided Design of Integrated Circuits and Systems

Keywords: Multicore architectures, task mapping, real-time and energy constraints, QoS, MILP, problem decomposition

Multicore architectures have great potential for energy-constrained embedded systems, such as energy-harvesting wireless sensor networks. Some embedded applications, especially the real-time ones, can be modeled as imprecise computation tasks. A task is divided into a mandatory subtask that provides a baseline Quality-of-Service (QoS) and an optional subtask that refines the result to increase the QoS. Combining dynamic voltage and frequency scaling, task allocation and task adjustment, we can maximize the system QoS under real-time and energy supply constraints. However, the nonlinear and combinatorial nature of this problem makes it difficult to solve. This work first formulates a mixed-integer non-linear programming problem to concurrently carry out task-to-processor allocation, frequencyto-task assignment and optional task adjustment. We provide a mixed-integer linear programming form of this formulation without performance degradation and we propose a novel decomposition algorithm to provide an optimal solution with reduced computation time compared to state-of-the-art optimal approaches (22.6% in average). We also propose a heuristic version that has negligible computation time.

I. INTRODUCTION

The increasing performance requirements of applications have pushed the embedded systems towards multicore architectures. Within thirty years, the code size of automotive, space and avionics applications has significantly increased [START_REF] Girbal | On the convergence of mainstream and mission-critical markets[END_REF]. Multicore architectures provide significant Space, Weight and Power savings (SWaP) while offering massive computing capabilities compared with single core processors. They are also capable of integrating diverse applications on the same platform [START_REF] Lemonnier | Towards future adaptive multiprocessor systems-on-chip: an innovative approach for flexible architectures[END_REF], [START_REF] Chen | Narrowband internet of things: Implementations and applications[END_REF].

In some application domains, e.g., audio/video streaming, radar tracking, and control system design, less accurate results computed before the deadline are preferable than accurate, but too late, results [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF]. This statement holds because a real-time application has to provide a result before a given deadline. When not enough time is available, approximate results are acceptable as long as the baseline Quality-of-Service (QoS) is satisfied and the results are provided in time [START_REF] Chishiro | Practical imprecise computation model: theory and practice[END_REF]. For instance, in audio/video streaming, frames with a lower quality are better than missing frames. In radar tracking, an estimation of target's location in time is better than an accurate location arriving too late. In control loops, an approximate result produced by a control law is more preferable as long as the controlled system, e.g., cruise control system, remains stable. In these domains, the applications can be modeled as L. Mo. A. Kritikakou, and O. Sentieys are with Univ Rennes, INRIA, IRISA, CNRS, 35042 Rennes Cedex, France. E-mail: lei.mo@inria.fr, angeliki.kritikakou@irisa.fr, olivier.sentieys@irisa.fr.

This article was presented in the International Conference on Hardware/Software Codesign and System Synthesis 2018 and appears as part of the ESWEEK-TCAD special issue.

Imprecise Computation (IC) tasks [START_REF] Liu | Imprecise computations[END_REF], where a task is logically decomposed into a mandatory subtask and an optional subtask. All the mandatory subtasks must be completed before the deadline to have an acceptable result, while the optional subtasks can be left incomplete at the cost of reduced quality. The QoS of such systems increases with the longer execution of the optional subtasks.

At the same time, the energy consumption has become an important concern, especially for the systems with limited energy budget, such as battery powered devices and energyharvesting systems. Dynamic Voltage and Frequency Scaling (DVFS) is a technique that can manage both voltage and frequency, and, thus, control the energy and the time required to execute the tasks. By adequately allocating the tasks onto the processors and deciding the voltage/frequency of the processor for each task execution, the QoS can be further improved under the same energy supply and real-time constraints.

The majority of the task mapping approaches on multicores, e.g., [START_REF] Mahmood | Energyaware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]- [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF], has been conducted on how to allocate the tasks and/or decide the voltage/frequency in order to reduce the energy consumption of the system and guarantee the real-time constraints, i.e., energy-aware task mapping. However, most of these works focus only on the execution of the mandatory tasks. Some approaches, e.g., [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise computation systems with energy considerations[END_REF]- [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF], aim to maximize QoS by exploring the execution of the optional subtasks, with energy supply and real-time constraints, i.e., QoS-aware task mapping. However, the aforementioned approaches focus on different contexts compared to the problem considered in this paper, as summarized in Table I. Besides single-objective optimization, other existing approaches focus on bi-objective optimization, e.g., minimizing the energy consumption and maximizing the QoS, but without restrictions on the energy supply [START_REF] Micheletto | Using bioinspired metaheuristics to solve reward-based energy-aware mandatory/optional realtime tasks scheduling[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. In addition, the task mapping problems on multicore platforms usually have an N P-hard complexity [START_REF] Fisher | The multiprocessor real-time scheduling of general task systems[END_REF]. Therefore, most of these QoS-aware task mapping approaches focus on heuristics to find feasible solutions [START_REF] Yu | Dynamic scheduling of imprecisecomputation tasks in maximizing QoS under energy constraints for embedded systems[END_REF]- [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF], [START_REF] Micheletto | Using bioinspired metaheuristics to solve reward-based energy-aware mandatory/optional realtime tasks scheduling[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF].

Complementary to the state-of-the-art, our work jointly optimizes task-to-processor allocation, frequency-to-task assignment and optional subtask adjustment, under energy supply and real-time constraints, with an objective of maximizing the QoS. Two algorithmic solutions are proposed: 1) an optimal approach, named Optimal Joint DVFS Task Mapping (OJTM), with lower computational complexity than the state-of-the-art optimal approaches, and 2) a polynomial-time complexity algorithm, named Heuristic Joint DVFS Task Mapping (HJTM), with better quality results compared with the state-of-theart heuristics. Heuristics are of great practical significance, since they are able to find feasible solutions in a short time. However, it is also important to find the optimal solution. Only by doing so, we can find out how far the obtained feasible solution is from the optimal one, and how to improve the heuristics based on the optimal solution. A. Related Work 1) Energy-aware Task Mapping: Existing works that focus on the pure energy-aware task mapping problem aim at minimizing the energy consumption under resource and application constraints, e.g., [START_REF] Mahmood | Energyaware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]- [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF]. DVFS is widely used in energy-aware task mapping. When voltage level is discrete, the problems are usually formulated as Integer Programming (IP), e.g., [START_REF] Mahmood | Energyaware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]- [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. To efficiently solve the IP problems, a hybrid Genetic Algorithm (GA) is proposed in [START_REF] Mahmood | Energyaware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF], a polynomialtime complexity two-step heuristic is designed in [START_REF] Xu | Energy minimizing for parallel realtime tasks based on level-packing[END_REF], and the IP problem is relaxed to a Linear Programming (LP) in [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. Combining DVFS and Dynamic Power Management (DPM), a Mixed-Integer Linear Programming (MILP)-based task mapping problem is considered in [START_REF] Chen | Energy optimization for realtime multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF] and the problem is solved using the CPLEX optimization solver. When the voltage level is continuous, a convex problem is formulated in [START_REF] Kong | Energy-efficient scheduling of real-time tasks on cluster-based multicores[END_REF] and the problem can be solved by using polynomial-time methods. In [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF], Mixed-Integer Non-Linear Programming (MINLP) is used to formulate the task mapping problem. The problem is relaxed to an MILP by linear approximation and solved by Branch and Bound (B&B) method [START_REF] Boyd | Branch and bound methods[END_REF]. However, the tasks in aforementioned approaches are not modifiable, and, thus, no exploration of the QoS improvement through the optional subtasks adjustment is considered.

2) QoS-aware Task Mapping: Other works consider the QoS-aware task mapping problem using the IC task model and having a goal to maximize the QoS under a set of realtime and/or energy supply constraints, e.g., [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise computation systems with energy considerations[END_REF]- [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. The target platforms considered in [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise computation systems with energy considerations[END_REF], [START_REF] Yu | Dynamic scheduling of imprecisecomputation tasks in maximizing QoS under energy constraints for embedded systems[END_REF] are single core platforms. Therefore, the task allocation problem is not taken into account. Although the works in [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Yu | Dynamic scheduling of imprecise computation tasks on real-time embedded multiprocessors[END_REF]- [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] target at multicore platform, some assumptions are made during the problem formulation. More precisely, the task-to-processor allocation is fixed and given in advance for all the tasks in [START_REF] Yu | Dynamic scheduling of imprecise computation tasks on real-time embedded multiprocessors[END_REF], while in [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF]- [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF] each processor has a predefined frequency, and in [START_REF] Micheletto | Using bioinspired metaheuristics to solve reward-based energy-aware mandatory/optional realtime tasks scheduling[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] the energy supply constraints are not taken into account since the aim is to maximize the QoS as well as to minimize the energy consumption. When taking multiple system requirements into account, the complex coupling between the optimization variables makes the problem difficult to solve, especially when the coupling is non-linear and non-convex. The methods that used to solve the aforementioned problems can be classified into two main classes. The first class includes the methods based on heuristics, e.g., [START_REF] Yu | Dynamic scheduling of imprecisecomputation tasks in maximizing QoS under energy constraints for embedded systems[END_REF]- [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF], [START_REF] Micheletto | Using bioinspired metaheuristics to solve reward-based energy-aware mandatory/optional realtime tasks scheduling[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. The second class includes the methods that always produce an optimal solution, e.g., [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise computation systems with energy considerations[END_REF], [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF].

B. Contributions

This paper solves the following problem: Given a DVFSenabled multicore platform and a set of IC tasks, we determine which processor should the task be executed on as well as the frequency of each task, such that the system QoS is maximized while guaranteeing the system constraints, i.e., energy supply and task deadline. Hence, we need to decide task-to-processor allocation, frequency-to-task assignment and task adjustment. Our main contributions are summarized as follows:

1) We propose an MINLP formulation of the QoS-aware mapping problem of IC tasks on DFVS multicore platforms under a set of energy supply and real-time constraints. The task-to-processor allocation, frequency-totask assignment, and task adjustment are optimized simultaneously. Solving these correlated subproblems separately may lead to non-optimal solutions. Then, applying the variable replacement method, the N Phard joint-design problem is safely transformed to an MILP. We prove that the proposed linearization does not cause performance degradation, compared with the existing linear approximation methods [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF]. 2) We design an OJTM solving approach based on Benders decomposition [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF]. The MILP-based task mapping problem is divided into two subproblems with simpler structures, i.e., with less constraints and variables. The first subproblem is an Integer Linear Programming (ILP) and it is responsible for the taskto-processor allocation and the frequency-to-task assignment. The second subproblem is an LP and it is related to the task adjustment. The solution is found through iterations. During one iteration, the solution of the first subproblem is used as input to the second sub-problem. The solution of the second subproblem is incorporated to the first one by adding new constraints for the next iteration. We prove that by iterating these two subproblems, the OJTM is guaranteed to reach the global optimal solution. 3) We present a novel HJTM algorithm that reduces the complexity of OJTM by removing the iteration process.

As the HJTM is a reduced version of OJTM, it is easy to be applied to similar MILP problems, compared with the existing heuristics. When the problem formulation is changed, existing heuristics usually have to be redeveloped, such as [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. 4) We evaluate the performance of the proposed approaches by a set of extended experiments, where the proposed approaches are compared both in QoS and in computation time with the state-of-the-art optimal, stochastic and heuristic approaches.

C. Paper Organization

The remainder of this paper is organized as follows: Section II presents the system model and the problem formulation. Section III describes the OJTM algorithm to optimally map IC tasks to DVFS-enabled multicore platform and in Section IV we present the details of HJTM algorithm. Section V evaluates the performance of OJTM and HJTM algorithms through several simulations and Section VI concludes this study.

II. MINLP PROBLEM FORMULATION

A. Task model

The applications of our target domain are modeled as IC tasks [START_REF] Liu | Imprecise computations[END_REF].

Definition 2.1 (Imprecise computation task): A task τ i is defined as an imprecise computation task when it can be logically divided into a mandatory subtask and an optional subtask. The mandatory subtask executes M i cycles to generate the baseline QoS level. The optional subtask executes o i cycles to improve the QoS.

Consider a task set of N independent IC tasks {τ 1 , . . . , τ N }. Task τ i is described by a tuple {o i , M i , O i , D i }, where O i is the upper bound on the number of possible optional cycles, i.e., 0 ≤ o i ≤ O i . Therefore, the total length of task τ i , measured in execution cycles, is M i + o i . The exact number of mandatory cycles M i and maximum optional cycles O i of a task τ i depends on the actual execution of the task on the multicore system [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Chishiro | Practical imprecise computation model: theory and practice[END_REF]. As we require to guarantee the deadline of the task, the M i and O i are measured in Worst Case Execution Cycles (WCECs) [START_REF] Yu | Dynamic scheduling of imprecisecomputation tasks in maximizing QoS under energy constraints for embedded systems[END_REF]. For each task τ i , there is a relative deadline D i before which both the mandatory and the optional subtasks of τ i must be completed. The relative deadline is defined as the time interval between the start of a task and the time instance when the deadline occurs [START_REF] Huang | Self-suspension real-time tasks under fixed-relative-deadline fixed-priority scheduling[END_REF]. The scheduling horizon H is given by the global period of the tasks [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF]. Hence, the relative deadline D i cannot be greater that the scheduling horizon H, i.e., D i ≤ H.

B. Energy model

We assume that the IC tasks run on a multicore platform that contains M homogeneous processors {θ 1 , . . . , θ M }. The processors can support DVFS and each processor

θ k has L differ- ent discrete voltage/frequency levels {(v 1 , f 1 ), . . . , (v L , f L )}.
We consider that the processor energy consumption is not only determined by the dynamic power, but also by the static power due to CMOS device feature scaling. The dynamic power by executing task τ i is related to processor clock frequency f l and supply voltage v d , expressed as

P d l = C ef f v 2 l f l ,
where C ef f is the average effective switching capacitance. For static power, we adopt P s l = v l K 1 e K2v l e K3v bs + |v bs |I j from [START_REF] Martin | Combined dynamic voltage scaling and adaptive body biasing for lower power microprocessors under dynamic workloads[END_REF]. The constants K 1 , K 2 and K 3 depend on the processor technology. I j is the approximately constant junction leakage current. v bs is a reverse bias voltage used to reduce the leakage power and can be treated as constant [START_REF] Jejurikar | Leakage aware dynamic voltage scaling for real-time embedded systems[END_REF]. Hence, the power of the processor under the voltage/frequency level (v l , f l ) is

P c l = P s l + P d l . ( 1
)
This model is widely adopted by the similarly published works [START_REF] Chen | Energy optimization for realtime multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], [START_REF] Yu | Dynamic scheduling of imprecise computation tasks on real-time embedded multiprocessors[END_REF], [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF].

For readability reasons in the rest of the manuscript, we use the term "frequency" instead of "voltage/frequency". When the frequency is determined, the corresponding voltage is also given. Examples of processor frequency-voltage characteristics are provided in Table II of the experimental section. In this paper, we consider inter-task DVFS [START_REF] Chen | Energy optimization for realtime multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], [START_REF] Yu | Dynamic scheduling of imprecisecomputation tasks in maximizing QoS under energy constraints for embedded systems[END_REF], i.e., the frequency of a processor stays constant for the entire execution of a task. Unlike the big.LITTLE platforms [START_REF] Chwa | Energy and feasibility optimal global scheduling framework on big.LITTLE platforms[END_REF] that support task migration, we consider that a task starts execution on a processor and ends its execution on the same processor [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF], [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. We assume that the system is energyconstrained in the sense that it has a fixed energy budget E s that cannot be replenished during the scheduling horizon H [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF], [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF]. Taking the available energy E s into account, the system operation can be divided into three states: 1) Low: the supplied energy E s is insufficient to execute all the mandatory cycles {M 1 , . . . , M N }, 2) High: the supplied energy E s is sufficient to execute all the mandatory and optional cycles {M 1 +O 1 , . . . , M N +O N }, and 3) Medium: all the mandatory cycles are ensured to finish, while not all the optional cycles can complete their executions. In this paper, we focus on medium state.

C. Problem Formulation

The objective is to maximize the QoS subject to a set of real-time and energy constrains. In this context, we determine 1) on which processor should the task be executed (task-toprocessor allocation), 2) what frequency should be used for the task (frequency-to-task assignment), and 3) how many optional cycles should be executed (task adjustment). Therefore, we define the following variables:

1) Binary variables

q ik = 1 if task τ i runs on processor θ k , 0 otherwise. c il = 1 if task τ i executes with frequency f l , 0 otherwise. 2) Continuous variable o i optional cycles of task τ i .
For tractability reasons, o i is considered as continuous variable in the problem formulation. When the problem is solved, the result is rounded down. As the tasks execute typically hundreds of thousands of cycles, this impact is negligible [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. Let N {1, . . . , N }, M {1, . . . , M } and L {1, . . . , L}. Since 1) each task can be assigned only one frequency level, and 2) each task is executed at only one processor, the binary variables c il and q ik must satisfy the following constraints:

L l=1 c il = 1, ∀i ∈ N , (2)
M k=1 q ik = 1, ∀i ∈ N . (3)
For the real-time constraints, since 1) the execution time of task M i + o i cannot exceed the relative deadline D i , and 2) the tasks assigned to processor θ k must be executed within the scheduling horizon H, we have

L l=1 c il w il + o i f l ≤ D i , ∀i ∈ N , (4)
N i=1 q ik L l=1 c il w il + o i f l ≤ H, ∀k ∈ M, (5)
where w il = Mi f l is the execution time of the mandatory subtask M i with frequency f l . D i is determined by the applications and

D i ≥ min ∀l∈L Mi f l .
Based on the energy model (1), the total energy consumed by M processors within the horizon H is bounded by

N i=1 L l=1 c il w il + o i f l P s l + P d l + M H - N i=1 L l=1 c il w il + o i f l P s 0 ≤ E s , (6)
where M H -

N i=1 L l=1 c il w il + oi f l
is the time of M processors in the idle state during the scheduling horizon H and P s 0 is the static power of idle state. It has been shown in [START_REF] Aydin | Optimal reward-based scheduling for periodic real-time tasks[END_REF], [START_REF] Yu | Dynamic scheduling of imprecise computation tasks on real-time embedded multiprocessors[END_REF], [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF], [START_REF] Yu | Quality-driven dynamic scheduling for real-time adaptive applications on multiprocessor systems[END_REF] that applications exist, e.g., Tier-2 coder in M-JPEG2000, where the QoS highly depends on its optional subtask and it can be represented as a linear function of the optional cycles. The linear function increases its value uniformly with the optional cycles. To quantify the relationship between QoS and optional cycles, each task τ i is associated with a linear function g i (o i ). In order to maximize QoS function N i=1 g i (o i ) (or to minimize its negative), based on all the aforementioned constraints, the Primal Problem (PP) is formulated as 3), ( 4), ( 5), [START_REF] Liu | Imprecise computations[END_REF],

PP : min c,q,o Φ = - N i=1 g i (o i ) s.t.    (2), (
c il , q ik ∈ {0, 1}, 0 ≤ o i ≤ O i , ∀i ∈ N , ∀l ∈ L, ∀k ∈ M.
Due to the products of optimization variables c il o i and q ik c il o i in the constraints (4), ( 5) and ( 6), the PP is an MINLP problem, which is difficult to directly solve because of the nonlinear items [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF].

Theorem 2.1: The MINLP-based QoS-aware IC-tasks mapping problem, i.e., the PP, is N P-hard.

Proof: Due to page limitations, the proof is excluded from the manuscript. However, the details can be found in [START_REF] Samuel | Non-convex mixed-integer nonlinear programming: a survey[END_REF].

III. OPTIMAL APPROACH

In this section, we present our optimal approach for solving the problem defined in Section II.

A. MILP formulation

We adopt the idea of variable replacement to efficiently solve the PP. In this way, the MINLP can be safely transformed to an MILP with a simpler structure, and, thus, easier to solve. To achieve this, we first introduce an intermediate (continuous) variable

t i = L l=1 c il w il + o i f l , ∀i ∈ N , (7)
which is the execution time of task τ i . Then, we linearize the nonlinear items c il o i and q ik t i according to Lemma 3.1. Lemma 3.1: Given constants s 1 > 0 and s 2 > 0 and two constraint spaces

S 1 = {[h, b, x]|h = bx, -s 1 ≤ x ≤ s 2 , b ∈ {0, 1}} and S 2 = {[h, b, x]|-bs 1 ≤ h ≤ bs 2 , h+bs 1 -x-s 1 ≤ 0, h -bs 2 -x + s 2 ≥ 0, b ∈ {0, 1}}, then S 1 S 2 . Proof: Based on h = bx and -s 1 ≤ x ≤ s 2 , we have -bs 1 ≤ h ≤ bs 2 . And further, we obtain (b -1)(x + s 1 ) ≤ 0 and (b -1)(x -s 2 ) ≥ 0 due to -s 1 ≤ x ≤ s 2 and b ∈ {0, 1}. Hence, we have h+bs 1 -x-s 1 ≤ 0 and h-bs 2 -x+s 2 ≥ 0. S 1 ⇒ S 2 holds. If b = 0, based on -bs 1 ≤ h ≤ bs 2 , h + bs 1 -x -s 1 ≤ 0 and h -bs 2 -x + s 2 ≥ 0, we have h = 0 and -s 1 ≤ x ≤ s 2 . Similarly, if b = 1, we have -s 1 ≤ h = x ≤ s 2 . S 2 ⇒ S 1 holds.
Although the lower bound of

t i is min ∀l∈L Mi f l ≥ 0 (= 0 when M i = 0), we assume that 0 ≤ t i ≤ D i since the range 0, min ∀l∈L Mi f l
is excluded by [START_REF] Mahmood | Energyaware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]. Note that c il ∈ {0, 1},

q ik ∈ {0, 1}, 0 ≤ o i ≤ O i and 0 ≤ t i ≤ D i .
We introduce two intermediate (continuous) variables h il and d ik to replace the nonlinear items c il o i and q ik t i , respectively. Therefore, the constraints ( 4) and ( 6) are safely replaced by

L l=1 c il w il + h il f l ≤ D i , ∀i ∈ N , (8)
N i=1 L l=1 c il w il + h il f l P d l + P s l -P s 0 + M HP s 0 ≤ E s , (9) h il ≤ c il O i , ∀i ∈ N , ∀l ∈ L, (10)
h il ≤ o i , ∀i ∈ N , ∀l ∈ L, (11)
o i -h il -(1 -c il ) O i ≤ 0, ∀i ∈ N , ∀l ∈ L. ( 12
)
Based on the intermediate variables t i and h il , the constraint (5) is safely replaced by

t i = L l=1 c il w il + h il f l , ∀i ∈ N , (13)
N i=1 d ik ≤ H, ∀k ∈ M, (14)
d ik ≤ q ik D i , ∀i ∈ N , ∀k ∈ M, (15)
d ik ≤ t i , ∀i ∈ N , ∀k ∈ M, (16)
t i -d ik -(1 -q ik ) D i ≤ 0, ∀i ∈ N , ∀k ∈ M, (17)
where ( 10)-( 12) and ( 15)-( 17) are the additional constraints introduced by the linearization. Consequently, the PP is rewritten as the MILP problem:

PP1 : min c,q,o, t,h,d Φ = - N i=1 g i (o i ) s.t.         
(2), ( 3), ( 8)

-(17), c il , q ik ∈ {0, 1}, 0 ≤ o i , h il ≤ O i , min ∀l∈L Mi f l ≤ t i ≤ D i , 0 ≤ d ik ≤ D i , ∀i ∈ N , ∀l ∈ L, ∀k ∈ M.
Remark 3.1: Lemma 3.1 implies that the variable replacement does not change the feasible region of the problem (S 1 S 2 ). Since the objective functions of PP and PP1 are the same, solving PP1 is equivalent to solving PP, i.e., the optimal objective function values of PP1 and PP are the same.

B. Optimal Joint DVFS Task Mapping

We propose an optimal algorithm, i.e., OJTM, to solve PP1 based on its special characteristics, i.e., the binary and continuous variables are coupled with each other linearly.

The key idea of OJTM comes from Benders decomposition. Instead of considering all the variables and the constraints of PP1 simultaneously, OJTM decomposes the PP1 into a Master Problem (MP) and a Slave Problem (SP), which are solved iteratively through a feedback loop. By doing so, the computational complexity of the solution is significantly reduced [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF], [START_REF] Mcdaniel | A modified Benders' partitioning algorithm for mixed integer programming[END_REF], even if the PP1 is non-convex. The structure of OJTM is shown in Fig. 1. The MP accounts for all the binary variables (task-to-processor allocation and frequency-to-task assignment) and the corresponding part of the objective function and the constraints of PP1. It also includes a set of constraints called Benders cuts. The SP includes all the continuous variables (task adjustment) and the associated constraints of PP1. The SP solution provides additional constraint, which is incorporated in the MP through the Benders cuts.

Task-to-processor allocation

Frequency-to-task assignment For presentation convenience, matrices and vectors are used to denote the constraints and variables. Hence, the PP1 is reformulated as

PP1 : min f y s.t. Ax b 1 , Cx + Dy b 2 ,
where x is a binary vector, and y is a continuous vector. Vector f represents the objective function coefficients, and f represents the transposition of f . Matrices A, C and D represent the constraints' coefficients. b 1 and b 2 are the udimensional and v-dimensional vectors, respectively. Based on the structure of the PP1, at the m th iteration, the MP is formulated as

MP : Φ l (m) = min x, Φ Φ (18)
s.t.

     Ax b 1 , Φ ≥ λ(ς) (Cx -b 2 ), ∀ς ∈ A, 0 ≥ λ(ϑ) (Cx -b 2 ), ∀ϑ ∈ B,
where where x(m) is the MP solution at the m th iteration. In fact, the SP is "identical" to PP1: the formulations are the same, except that the binary variables x(m) in the SP are fixed. Denote Φ * as the optimal objective function value of the PP1. Note that the MP [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF] only contains the frequency-to-task and task-to-processor information, i.e., x(m). Compared with the PP1, the task adjustment related constraints are relaxed. Solving the MP yields a lower bound of Φ * , i.e., Φ l (m). On the other hand, x(m) may be just a feasible solution (not optimal yet), and SP is a minimization problem. Solving the SP with x(m) yields an upper bound of Φ * , i.e., Φ u (m) (the proof is provided in Theorem 3.1). Therefore, we have Φ l (m) ≤ Φ * ≤ Φ u (m). At each iteration, introducing a new feasibility constraint [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] or infeasibility constraint [START_REF] Boyd | Branch and bound methods[END_REF] into the MP, the gap between the upper and lower bounds, i.e., Φ u (m) -Φ l (m), gradually reduces (the proof is provided in Lemma 3.3), and the iterations stop when Φ u (m)-Φ l (m) ≤ ε.

Φ ≥ λ(ς) (Cx -b 2 ) (∀ς ∈ A) and 0 ≥ λ(ϑ) (Cx - b 2 ) (∀ϑ ∈ B)
Definition 3.1 (ε-optimal solution): A solution ( x, ỹ) of the PP1 is an ε-optimal solution, if it satisfies the constraint 0 ≤ f ỹ -Φ * ≤ ε.

Lemma 3.2: Based on the iteration termination criterion Φ u (m) -Φ l (m) ≤ ε, the solution (x(m), y(m)) obtained by OJTM is an ε-optimal solution.

Proof: Substituting the solution (x(m), y(m)) into the PP1, we have

Φ l (m) ≤ f y(m) = Φ u (m). Since Φ u (m) - Φ l (m) ≤ ε, we get 0 ≤ f y(m) -Φ * ≤ ε, which means (x(m), y(m)) is an ε-optimal solution.
If ε is a small positive value, the optimal value Φ * is found by OJTM (the proof is provided in Theorem 3.2). The iteration process is summarized as follows:

1) Step 1 -Initialization: Initialize the iteration counter m = 0, the solution x(0) of the MP, the lower bound Φ l (0) = -∞, and the upper bound Φ u (0) = +∞. The feasibility and infeasibility constraints are set to null. Note that the initial solution x(0) can be given arbitrarily, as long as it satisfies the constraints (2) and (3).

2) Step 2 -Solving SP: Assume that the current iteration is m. Since the SP is a LP problem, the optimal objective function values of the SP and its dual problem are equivalent due to the strong duality [START_REF] Boyd | Convex optimization[END_REF]. Introducing the positive Lagrange multipliers λ [λ i ] (1 ≤ i ≤ v) to the SP, the associated dual problem (DSP) is

DSP : max λ 0 λ (Cx(m) -b 2 ) (19) s.t. f + D λ 0.
the DSP is a LP problem, the solution can be found quickly using standard algorithms, such as simplex method or interior point method [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF].

3) Step 3 -Solving MP: Depends on the following cases: 1) The DSP is infeasible: The SP has an unbounded solution. Hence, the PP1 has no feasible solution.

2) The DSP has a bounded solution λ(m): The SP is feasible due to the strong duality, and A ← {m} ∪ A.

Assume that Φ(m) and y(m) are the solution of the MP and the SP at the m th iteration, respectively. Since i) x(m) is a feasible solution (not optimal) of PP1, and ii) f y(m) = λ(m) (Cx(m) -b 2 ) due to the strong duality, the upper bound of Φ * at the m th iteration is updated by

Φ u (m) = min{Φ u (m - 1), λ(m) (Cx(m) -b 2 )}. Note that i) Φ(m) < λ(m) (Cx(m) -b 2 )
, and ii) x(m) is not the optimal solution. In order to avoid selecting the non-optimal solution x(m) again, a new feasibility constraint

Φ ≥ λ(m) (Cx -b 2 ) (20)
is generated and added into the MP at the (m + 1) th iteration.

3) The DSP has an unbounded solution, i.e., λ(m) (Cx(m) -b 2 ) = +∞: Due to the strong duality, the SP has no feasible solution under given solution x(m), and B ← {m}∪B. For the SP, its feasibility is related to the constraints rather than the objective function. This problem may feasible if the positive variables ξ [ξ i ] (1 ≤ i ≤ v) are introduced to relax the constraints.

Based on this idea, we construct a Feasibility Check Problem (FCP)

FCP : min ξ 0 1 ξ s.t. Cx(m) + Dy b 2 + ξ.
Since the FCP is a LP problem, the strong duality exists between the FCP and its dual problem. Introducing the positive Lagrange multipliers λ = λi (1 ≤ i ≤ v) to the FCP, the associated dual problem (DFCP) is

DFCP : max λ 0 λ (Cx(m) -b 2 ) (21)
s.t. 1 -λ 0.

The presented method used to solve DSP can be also applied here to solve FCP and DFCP, as both are LP problems. Assume that ξ(m) and λ(m) are the solutions of FCP and DFCP, respectively. If the SP exists infeasible constraints, the related relax variables are non-zero, while the others are zero, i.e., 1 ξ(m) > 0. Due to the strong duality, we have

1 ξ(m) = λ(m) (Cx(m) -b 2 ) > 0.
In order to avoid selecting the infeasible solution x(m) again, a new infeasibility constraint

0 ≥ λ(m) (Cx -b 2 ) (22)
is generated and added into the MP at the (m + 1) th iteration. Note that in [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] and ( 22), all the parameters are constant except Φ and x. These are the variables of the MP at the (m+ 1) th iteration. With this new added feasibility/infeasibility constraint, the MP is solved again to obtain a new solution x(m + 1) at the next iteration.

Remark 3.2: Note that 1) the SP has the same objective function as the PP1, i.e., f y, and 2) the SP and the DSP are equivalent due to the strong duality, i.e., they have the same optimal objective function values. From [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF], we can see that the auxiliary variable Φ has the same physical meaning as the objective function of PP1.

Remark 3.3: We solve the DSP rather than solving the SP directly, because 1) the SP and its dual problem are equivalent due to the strong duality, and 2) we can construct feasibility and infeasibility constraints through the solution of the DSP. The reason why we select the solution of the DFCB rather than the solution of the FCB to construct infeasibility constraint is that λ (Cx -b 2 ) is a function with respect to x but not 1 ξ, i.e., 0 ≥ 1 ξ is an invalid constraint for the MP.

C. Convergence analysis

Equation [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF] shows that 1) the MP is composed of the binary variables x and a continuous variable Φ, 2) the real constraints are Ax b 1 and 0 ≥ λ(ϑ) (Cx -b 2 ) (∀ϑ ∈ B), and 3) λ(ς) (Cx -b 2 ) (∀ς ∈ A) can be treated as objective functions. Therefore, the MP can be solved by only considering the binary variables x [START_REF] Mcdaniel | A modified Benders' partitioning algorithm for mixed integer programming[END_REF]. Comparing the MP with the following optimization problem, which has the same constraints, we have Φ(m) = max ∀i∈A Φr (i) .

MP i : Φr (i) = min x λ(i) (Cx -b 2 ) s.t. Ax b 1 , 0 ≥ λ(ϑ) (Cx -b 2 ), ∀ϑ ∈ B.
Theorem 3.1: Solving the MP and the SP at the m th iteration, we obtain a lower bound Φ l (m) and an upper bound Φ u (m) of the optimal value of the objective function Φ * , respectively.

Proof: First, we prove that Φ l (m) = Φ(m) is a lower bound of Φ * , where Φ(m) is the solution of the MP at the m th iteration. Without loss of generality, we assume that Φ(m) = Φr (ρ) = max ∀i∈A Φr (i) , where ρ ∈ A. Hence, we have

Φ(m) = min x λ(ρ) (Cx -b 2 ) ≤ λ(ρ) (Cx * -b 2 ) (23a) ≤ max λ 0 λ (Cx * -b 2 ) = Φ * , (23b)
where min x λ(ρ) (Cx-b 2 ) is the objective function value of the MP ρ , and max λ 0 λ (Cx * -b 2 ) is the objective function value of the DSP. (23a) holds since x * is not the optimal solution of problem min x λ(ρ) (Cx -b 2 ). (23b) holds since solving the DSP with the optimal binary variables x * we find optimal value Φ * . [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] shows that Φ l (m) = Φ(m) is a lower bound of Φ * .

In the following, we prove that

Φ u (m) = min{Φ u (m - 1), λ(m) (Cx(m) -b 2 )} is an upper bound of Φ * . Note that we have min{Φ u (m -1), λ(m) (Cx(m) -b 2 )} = min 1≤i≤m {λ(i) (Cx(i) -b 2 )},
where x(i) and λ(i) are the solutions of the MP and the DSP at the i th iteration, respectively. Depending on the solution of the DSP, the objective function value of the DSP, i.e., λ(i) (Cx(i) -b 2 ), can be either finite or infinite. If the DSP has an unbounded solution, i.e., λ(i) (Cx(i) -b 2 ) = +∞, it is obvious that +∞ is an upper bound of Φ * . Thus, we consider the case when the DSP has a bounded solution λ(i). Due to the strong SP and DSP duality, we have

λ(i) (Cx(i) -b 2 ) = min y 0 f y|x(i) ≥ min y 0 f y|x * = Φ * , 1 ≤ i ≤ m,
where min y 0 f y|x(i) is the objective function value of the DSP with x(i). Thus, Φ u (m) = min 1≤i≤m {λ(i) (Cx(i)b 2 )} is an upper bound of Φ * .

Lemma 3.3: The lower bound sequence {Φ l (m)} is increasing, while the upper bound sequence {Φ u (m)} is decreasing.

Proof: Assume that ( Φ(m), x(m)) is the solution of the MP at the m th iteration. If x(m) are the non-optimal (infeasible) values of variables x, x(m) are excluded by adding new constraint

Φ ≥ λ(m) (Cx -b 2 ) (0 ≥ λ(m) (Cx -b 2 )
) into the MP at the (m + 1) th iteration. Note that the variables of the MP, i.e., x and Φ, are coupled with each other linearly. If x(m) are excluded, Φ(m) is excluded as well. Since 1) the aim of the MP is to minimize the objective function, 2) the non-optimal values of Φ * , i.e., { Φ(0), . . . , Φ(m)}, have been excluded by the constraints and3) with iteration number m increasing, more constrains are added into the MP (i.e., the feasible region of the MP will shrink), Φ l (m+1) = Φ(m+1) is larger than the previous lower bounds {Φ l (0), . . . , Φ l (m)}.

Φ ≥ λ(ς) (Cx-b 2 ) (∀ς ∈ A) and 0 ≥ λ(ϑ) (Cx -b 2 ) (∀ϑ ∈ B),
Assume that Φ u (m + 1) > Φ u (m). This contradicts the fact that

Φ u (m + 1) ≤ Φ u (m) and Φ u (m + 1) ≤ λ(m + 1) (Cx(m+1)-b 2 ) due to Φ u (m+1) = min{Φ u (m), λ(m+ 1) (Cx(m+1)-b 2 )}. Therefore, Φ u (m+1) is no larger than the previous upper bounds {Φ u (0), . . . , Φ u (m)}.
Theorem 3.2: each iteration with a new feasibility constraint [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] or infeasibility constraint [START_REF] Boyd | Branch and bound methods[END_REF] added into the MP, the solution obtained by OJTM converges to the global optimal value within a finite number of iterations.

Proof: At each iteration, there is a new constraint ((20) or ( 22)) added into the MP to exclude the non-optimal or the infeasible values of the binary variables x. Note that 1) the lower bound sequence {Φ l (0), . . . , Φ l (m)} is increasing, 2) the upper bound sequence {Φ u (0), . . . , Φ u (m)} is decreasing, and 3) the dimension of binary variables x is finite. According to the Theorem 2.4. in [START_REF] Geoffrion | Generalized Benders decomposition[END_REF], the solution converges to the global optimal value within a finite number of iterations.

IV. HEURISTIC APPROACH

The computational complexity of OJTM is dominated by the cost of solving the ILP-based MP at each iteration. Although the solution provided by OJTM is optimal, this method cannot be used to efficiently solve large problem sizes. To enhance the scalability of the proposed approach, we provide a novel heuristic algorithm, i.e., HJTM, to efficiently solve the PP. The reasons why we target PP rather than PP1 are that 1) these two problems are equivalent, and 2) the PP contains less variables and constraints making the heuristic less complex. The main idea of the HJTM is that the most important step of solving the PP is to find proper task-to-processor and frequency-to-task decisions (similar to basic idea of the OJTM). If the binary variables c and q are determined, the PP reduces to a LP problem. The structure of HJTM is shown in Fig. 2. HJTM has three steps. During the first step, we obtain feasible frequency-to-task assignment. Having determined the frequency-to-task decision, the problem reduces to finding feasible task allocation to the processors. In the final step, the optional cycles are defined to meet the real-time and the energy supply constraints.

Task-to-processor decision

Frequency-to-task decision FAP (ILP) TAP (ILP)

OTAP (LP)

Task adjustment decision

Feasibility pump

Relaxed problem (LP) Fig. 2. The structure of HJTM algorithm.

A. Frequency-to-task Assignment

Since the mandatory subtasks must be always executed, we initially consider only the frequency assignment of the mandatory subtasks and our aim is to minimize the energy consumption by assigning proper frequency to each task. From the constraint (6), we observe that if o i = 0 (∀i ∈ N ), the consumed energy of M processors reduces to

N i=1 L l=1 c il w il P s l + P d l -P s 0 + M HP s 0 . Since item M HP s 0 is fixed, the Frequency Assignment Problem (FAP) is formulated as FAP : min c N i=1 L l=1 c il w il P s l + P d l -P s 0 s.t.
(2),

L l=1 c il w il ≤ D i , ∀i ∈ N .

B. Task-to-processor Allocation

Solving the FAP, we obtain the execution time of task τ i 's mandatory subtask, i.e., L l=1 c il Mi f l . Taking task allocation decision q into account, the execution time of processor θ k is computed by

Λ(k) = N i=1 q ik L l=1 c il Mi f l .
Since the processors have the same voltage/frequency levels {(v 1 , f 1 ), . . . , (v L , f L )}, different task allocation schemes will not influence the solution of FAP, i.e., the task-to-processor decision can be made when frequency-to-task assignment is done. We introduce an auxiliary(continuous) variable U , which is denoted as the maximum task execution time among the processors, i.e., U = max ∀k∈M {Λ(k)}, and our aim is to minimize U by allocating the selected tasks to each processor. Thus, the Task Allocation Problem (TAP) is formulated as

TAP : min q,U U s.t.          U ≥ N i=1 q ik L l=1 c il M i f l , ∀k ∈ M, (3),
N i=1 q ik L l=1 c il M i f l ≤ H, ∀k ∈ M.
Note that the TAP is composed of binary variables q and a continuous variable U . It can be solved using a similar method as the one presented for the solution of the MP.

C. Optional Task Adjustment

Based on frequency-to-task assignment decision c and taskto-processor allocation decision q, the PP reduces to the following Optional Task Adjustment Problem (OTAP): 4), ( 5), [START_REF] Liu | Imprecise computations[END_REF] under given c and q.

OTAP : min o - N i=1 g i (o i ) s.t. (
Although the FAP and the TAP are ILP, compared with LP, they are still complex to solve directly, especially when problem size is large. To reduce the computational complexity, we relax FAP and TAP to LP problems. In particular, for the FAP (TAP), we introduce continuous variables q (ĉ), where 0 ≤ qik ≤ 1 (0 ≤ ĉil ≤ 1), to replace the original integer variables q (c). And then, we solve the relaxed problem (LP) and round the solution q (ĉ) to the nearest binary matrix that is feasible to the FAP (TAP). Such a binary matrix can be efficiently found by using the Feasibility Pump (FP) method proposed in [START_REF] Matteo | The feasibility pump[END_REF].

V. EXPERIMENTAL RESULTS

The values and the tuned parameters of the experimental set-up are summarized in Table II. The processor parameters are adopted from [START_REF] Chen | Energy optimization for realtime multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], where the platform is based on 70 nm technology and the accuracy of the parameters' values has been verified by SPICE simulations. The number of processors M is tuned from 4 to 10 with a step of 2. The WCECs of tasks are assumed be in the range provided in Table II. This range is calculated in [START_REF] Ramkumar | Temperature aware task sequencing and voltage scaling[END_REF] from the MiBench and MediaBench benchmark suites. The number of tasks N is tuned from 10 to 50 with a step of 10. The QoS function of the tasks, i.e., the objective function of the PP, is adopted from [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF], [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF]. The scheduling horizon H is assumed to be proportional to the average processor workload, and the energy supply E s is set to E s = ηE h , where E h is the minimum energy required to execute {M 1 + O 1 , . . . , M N + O N } cycles, η ∈ [0, 1] is an energy efficiency factor, and η is tuned from 0.8 to 0.9 with a step of 0.05. Note that different processor and task parameters lead to different values in the parameters {A, C, D, f , b 1 , b 2 } for the PP1. However, the problem structures under different values of parameters are the same, and, thus, the proposed methods are still applicative.

M i +O i f l Objective function N i=1 gi(oi) = N i=1 oi Constraints H = N M N i=1 D i N Es = ηE h E h = M HP s 0 + N i=1 min ∀l∈L M i +O i f l P s l + P d l -P s 0 Tuned parameters M N η Min/Max/Step 4/10/2
10/50/10 0.8/0.9/0.05

We compare the behavior (QoS and computation time) of the proposed OJTM and HJTM with: 1) optimal approaches, i.e. Branch and Bound method (B&B) [START_REF] Boyd | Branch and bound methods[END_REF], which is known to provide optimal solution for the MILP problem -as far as we know no optimal algorithm exist for the problem formulation PP1, 2) stochastic approaches, i.e. Genetic Algorithm (GA) [START_REF] Sayuti | Real-time low-power task mapping in networks-on-chip[END_REF], and 3) heuristics, i.e. Adaptive Task Allocation (ATA) [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF] -a two-step optimization approach. The simulations are performed on a laptop with quad-core 2.5 GHz Intel i7 processor and 16 GB RAM, and the algorithms are implemented in Matlab 2016a.

The QoS achieved by OJTM for all M , N and η parameters is shown in Fig. 3. During the simulations, the task sets in the experiments with N = i (i = 20, . . . , 50) and N = i -1 are correlated, i.e., the task set in the experiment with N = i is extended based on the task set in the experiment with N = i -1. Note that 1) the scheduling horizon H is proportional to the average processor workload, 2) the energy supply E s is proportional to the number of tasks N , and 3) the QoS is given by N i=1 o i . The QoS generally increases with respect to 1) N and 2) η under given M and N . We further explore the behavior of OJTM and HJTM by comparing the QoS gain of each method with the optimal B&B for all the aforementioned cases. Fig. 4 shows the statistical property of the QoS gain. The QoS gain of OJTM and HJTM (OJTM and B&B) is given by Qo(M,N,η)-Q h (M,N,η)

Qo(M,N,η) Qo(M,N,η)-Q b (M,N,η) Qo(M,N,η)
, where Q b (M, N, η), Q o (M, N, η) and Q h (M, N, η) are the QoS achieved by B&B, OJTM and HJTM under given M , N and η parameters, respectively. The box plot "OJTM vs HJTM" ("OJTM vs B&B") with η = 0.8 shows the statistical property of data set Qo(M,N,0.8)-Q h (M,N,0.8)

Qo(M,N,0.8) Qo(M,N,0.8)-Q b (M,N,0.8) Qo(M,N,0.8)
for all M and N parameters. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25 th and 75 th percentiles, respectively. The whiskers extend to the most extreme data points that are not considered outliers and the outliers are plotted individually using the '+' symbol. From Fig. 4, we observe that 1) the solutions given by B&B and OJTM are same, and, thus, OJTM also finds the optimal solution, and 2) OJTM achieves higher QoS (26.3% in average) than that of HJTM. The computation time of OJTM, HJTM and B&B is evaluated in Fig. 5. From Fig. 5(a), which shows the OJTM (HJTM) computation time for solving the PP1 (PP), we observe that 1) the HJTM computation time is almost unchanged with M , N and η varying, and 2) when M and N are fixed and η is tuned from 0.8 to 0.9, the changes of the OJTM computation time are small, compared with the changes of the OJTM computation time with M and N varying. Therefore, compared with OJTM, HJTM has a negligible computation time. We further compare the computation time gain of OJTM and B&B for all M , N and η parameters. The result is shown in Fig. 5 for all M and N parameters. With N increasing, the computation time of B&B and OJTM grows. However, OJTM takes a shorter computation time than that of B&B (22.6% in average). For an optimization problem, its computational complexity increases significantly with the number of variables and constraints. Solving iteratively smaller problems, i.e., the MP and the SP, is more efficient than solving a single large problem. This result agrees with the comparison of [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF]: the decomposition-based method is faster than the B&B for lager problem instances.

We define the OJTM convergence iteration as the number of iterations to achieve Φ u (m) -Φ l (m) ≤ ε, where ε = 0.01. Fig. 6 shows that with M and N increasing, more constraints are involved into the problem, and, thus, a higher number of iterations is required to converge to the optimal solution. We observe that when M ≤ 8 and N ≤ 40, the OJTM convergence iteration almost increases linearly with M and N .

Denote Q g (M, N, η) and T g (M, N, η) as the QoS achieved by GA and the computation time of GA under given M , N and η parameters, respectively. The QoS and computation time gains of OJTM and GA, i.e., , for all M , N η parameters, are shown in Fig. 7. From this figure we observe that 1) OJTM achieves a higher QoS than that of GA (7.6% in average), and 2) OJTM takes a shorter computation time than that of GA (35.6% in average). Although GA can solve complex nonlinear programming problem (non-convex), such as MINLP, the optimality of the solution is hard to guarantee. Compared with OJTM, the GA structure is more complex, as in each iteration GA needs to generate new populations. Therefore, the problem transformation from MINLP-based PP to MILPbased PP1 is necessary, since it can simplify the structure of the problem, and, thus, the optimal solution is much easier to find. Denote Q a (M, N, η) and T a (M, N, η) as the QoS achieved by ATA and the computation time of ATA under given M , N and η parameters, respectively. The QoS gain of HJTM and ATA, i.e., Q h (M,N,η)-Qa(M,N,η)

Q h (M,N,η)
, and the computation time gap of HJTM and ATA, i.e., {T h (M, N, η) -T a (M, N, η)}, for all M , N and η parameters, are shown in Fig. 8. ATA solves the task-to-processor allocation and the optional subtask adjustment problems in sequence, where each processor its own frequency. Since the HJTM optimizes also the frequencyto-task assignment, and, thus, achieves a higher QoS (96.7% in average) than that of ATA. We observe that 1) the difference between T h (M, η) and T a (M, N, η) is small, and 2) the average value of {T h (M, N, η)} is slightly higher than {T a (M, N, η)}. HJTM and ATA are both polynomial-time algorithms. Sometimes T h (M, N, η) can be a bit larger than T a (M, N, η), since 1) HJTM is a three-step heuristic while ATA is a two-step heuristic, and 2) we need to run FP method repeatedly until feasible solutions are found when solving the FAP and the TAP.

VI. CONCLUSION

In this paper, we formulated the problem of task mapping on DVFS-enabled homogeneous multicore as an MILP, with the goal of maximizing QoS without violating the real-time and the energy supply constraints. This problem is optimally solved by the OJTM algorithm. A novel algorithm, HJTM, is proposed to further reduce the computational complexity. Our numerical results show that OJTM is guaranteed to converge to the optimal solution, while HJTM can find a feasible solution within a negligible computation time, compared with OJTM. In the future, we plan to extend our method to heterogeneous multicore platforms (e.g., big.LITTLE platforms), accelerate the OJTM and consider the dependency among the tasks.

Fig. 1 .

 1 Fig. 1. The structure of OJTM algorithm.

  are the sets of feasibility and infeasibility constraints (Benders cuts), respectively. λ(m) and λ(m) are the solutions of the Dual Slave Problem (DSP)[START_REF] Micheletto | Using bioinspired metaheuristics to solve reward-based energy-aware mandatory/optional realtime tasks scheduling[END_REF] and the Dual Feasibility Check Problem (DFCP)[START_REF] Fisher | The multiprocessor real-time scheduling of general task systems[END_REF] at the m th iteration, respectively. A and B are the sets of iterations that the DSP has bounded and unbounded solutions, respectively. Note that 1) the objective function of the PP1 only contains continuous variables y, and 2) the MP only considers binary variables x. We introduce an auxiliary (continuous) variable Φ for the MP as the objective function, where Φ has the same physical meaning as the objective function of the PP1. The relationship between Φ and Φ is explained in the next section.The corresponding SP is formulated as SP : Φ u (m) = min y 0 f y s.t. Cx(m) + Dy b 2 ,

Fig. 3 .

 3 Fig. 3. QoS of OJTM with M , N and η varying.

Fig. 4 .

 4 Fig. 4. QoS gain of OJTM, HJTM and B&B with M , N and η varying.

  (b). Denote T b (M, N, η) and T o (M, N, η) as the computation time of B&B and OJTM under given M , N and η parameters, respectively. The box plot "B&B vs OJTM" with η = 0.8 the statistical property of data set T b (M,N,0.8)-To(M,N,0.8) T b (M,N,0.8)

  Qo(M,N,η)-Qg(M,N,η) Qo(M,N,η) and Tg(M,N,η)-To(M,N,η) Tg(M,N,η)

  Computation time of OJTM and HJTM with M , N and η varying.

  Computation time gain of OJTM and B&B with M , N and η varying.

Fig. 5 .

 5 Fig. 5. Comparison of OJTM, HJTM and B&B computation time with M , N and η varying.

Fig. 6 .

 6 Fig. 6. OJTM convergence iteration.

Fig. 7 .

 7 Fig. 7. QoS and computation time gains of OJTM and GA with M , N and η varying.

Fig. 8 .

 8 Fig. 8. QoS and computation time comparison of HJTM and ATA with M , N and η varying.

TABLE I

 I

					.		TASK MAPPING METHOD
			Energy-aware								QoS-aware
	[10]	[11]	[12]	[9]	[7]	[8]	[4]	[13]	[20]	[19]	[18]	[17]	[14]	[15]	[16]

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCAD.2018.2857300, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems

	3
	the-art heuristics, HJTM achieves an average 96.7%
	QoS improvement.
	The obtained re-
	sults show an average 22.6% reduction in the compu-
	tation time for the OJTM compared with state-of-the-
	art optimal approaches. On the other hand, HJTM runs
	∼100 times faster than OJTM with a cost of 26.3%
	QoS reduction. Moreover, compared with the state-of-

0278-0070 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

TABLE II

 II

		.	EXPERIMENTAL SET-UP
		Processor θ k characteristics
	v l (V)	0.65	0.7	0.75	0.8	0.85
	f l (GHz)	1.01	1.26	1.53	1.81	2.10
	P d l (mW) P s l (mW) P s 0 (µW)	184.9 246	266.7 370.4 498.9 290.1 340.3 397.6 80	655.5 462.7
		Task τi characteristics
	Mi, Oi ∈ [4 × 10 7 , 6 × 10 8 ]	Di = min ∀l∈L

√ √ √ √

ACKNOWLEDGMENT This research is funded by INRIA post-doctoral research fellowship program, and is partially funded by ANR ARTE-FACT (AppRoximaTivE Flexible Circuits and Computing for IoT) project (Grant No. ANR-15-CE25-0015), and is partly sponsored by the National Natural Science foundation of China (Grant No. 61403340).


Mohamed Raii

Francisco Javier

Escudero Sanz

Ange Nzihou

Rheological Behavior of Gypsum, Plaster, and Hydroxyapatite Gel Blends

The rheological behavior of water-calcium sulfates and hydroxyapatite-calcium sulfates binary blends was considered in this study. These blends are to be used as permeable reactive barriers for groundwater remediation. The influence of gypsum and plaster on the rheological behavior of Ca-HA Gel has been studied. The results show that all blends and formulations exhibit a shear-thinning effect and thixotropic behavior. The viscosity also depends on the level of concentration of the solid and the reactivity of the blends. The hydration of plaster causes cohesion and agglomeration between particles due to electrostatic forces, while the hydration of gypsum leads to dispersion due to repulsive forces between particles. The ζ potential was used in this study to understand the interaction between particles and its effect on the global behavior of the blends. The dissolution in water of sulfate ions containing calcium sulfate leads to a change in the surface charge of the blend. This allowed the dispersion or agglomeration of the mixtures to be predicted.

INTRODUCTION

Calcium phosphate is a well-known biomaterial frequently used in the medical domain and as an agricultural fertilizer. It can come from natural sources such as fish or cow bones or as a result of a synthesis process using calcium and phosphate compounds. Calcium phosphate has been investigated as a reactive material, especially in its thermodynamically stable structure with the stoichiometric ratio Ca/P = 1.67 (hydroxyapatite, Ca-HA), to determine its usefulness in treating contaminated water or as a stabilizing agent for hazardous waste. [START_REF] Chen | Evaluation of heavy metal remediation using mineral apatite[END_REF][START_REF] Nzihou | Calcium phosphate stabilization of fly ash with chloride extraction[END_REF][START_REF] Baillez | Removal of aqueous lead ions by hydroxyapatites: Equilibria and kinetic processes[END_REF] Many studies have also shown the capacity of Ca-HA to remediate soil and groundwater. [START_REF] Moore | Stabilization and Solidification of Hazardours[END_REF][START_REF] Conca | An Apatite II permeable reactive barrier to remediate groundwater containing Zn, Pb and Cd[END_REF] Gypsum can have two origins: natural gypsum and synthetic gypsum. Natural gypsum, characterized by a high purity, is an anionic mineral with a C 2h 6 space group. This gypsum is a rock resulting from the evaporation of water from saturated calcium sulfate and geological sedimentation. Synthetic gypsum is a byproduct of industrial activities. It is obtained from a reaction between sulfuric acid and phosphate rock to produce phosphoric acid (H 3 PO 4 ). The precipitate, usually named phosphogypsum (PG), is dihydrated calcium sulfate. This gypsum contains a significant amount of impurities and is considered to be a waste product. The environmental impact associated with PG storage and disposal, as reviewed by Tayibi et al., is negative on the surrounding land, water, and air. [START_REF] Tayibi | Environmental impact and management of phosphogypsum[END_REF] Valorization seems to be by far the best way to reuse this hazardous material. Industrial PG can be valorized as a material for soil stabilization. It can be added to Portland cement or to reactive materials such as fly ash to enhance the mechanical properties of the soil and to reduce the plasticity index, moisture, and water content. [START_REF] Degirmenci | Application of phosphogypsum in soil stabilization[END_REF] It can also be used as a hydraulic binder in the design of calcium sulfoaluminate cement. [START_REF] Kuryatnyk | Valorization of phosphogypsum as hydraulic binder[END_REF] The calcium sulfates used in this work are industrial gypsum (CaSO 4 •2H 2 O) and plaster (CaSO 4 •0.5H 2 O). They are byproducts of phosphoric acid production and are washed to eliminate the residual fraction (phosphate, fluor, sodium, radionuclides, etc.).

There are three anhydrous or anhydrite forms produced by sintering of hemihydrates or gypsum: the compact form α-CaSO 4 or CaSO 4 -I, obtained at 1200 °C and characterized by a crystalline structure; the stable form β-CaSO 4 or CaSO 4 -II, obtained at 400 °C, which is a natural insoluble anhydrite; the metastable form γ-CaSO 4 or CaSO 4 -III, obtained at low temperature between 130 and 200 °C, which is a soluble anhydrite. [START_REF] Iribarne | Reactivity of calcium sulfate from FBC systems[END_REF] The stability of gypsum is also influenced by the relative humidity (RH). At 10-35% RH, the migration of water molecules takes place toward the crystals, while at 100% RH, the migration takes place in the opposite direction. [START_REF] Auvray | The influence of relative humidity on the rate of convergence in an underground gypsum mine[END_REF][START_REF] Finot | Reactivity of gypsum faces according to the relative humidity by scanning force microscopy[END_REF] The crystalline structure of the gypsum, characterized by a perfect cleavage orientation (010), results in the circulation of water, leading to dissolution on the surface (010). [START_REF] Finot | Reactivity of gypsum faces according to the relative humidity by scanning force microscopy[END_REF][START_REF] Fan | Surface behavior of gypsum during dissolution[END_REF] The blending of hemihydrates with water gives a rigid and compact solid body, due to the hydration process.

The main purpose of this work is to investigate the rheological behavior of water/gypsum, water/plaster, Ca-HA/ gypsum, and Ca-HA/plaster blends. This is a preliminary step in a project devoted to the formulation and characterization of a Ca-HA-based permeable reactive barrier. The evaluation of the rheological behavior of a mixture of gypsum and plaster on Ca-HA Gel aims to better understand the development of the matrix microstructure and, consequently, to assess its workability. O, prepared by the wet method) used in this work are industrial byproducts of a phosphoric acid (H 3 PO 4 ) preparation from phosphate rock (PRAYON, Belgium). The plaster was used as received, and gypsum was sieved to under 1 mm and the larger fraction discarded.

MATERIALS AND METHODS

2.2. Hydroxyapatite Gel (Ca-HA Gel ). Ca-HA Gel was obtained by the reaction between calcite (CaCO 3 ; 98% pure, Fisher Scientific) and ammonium phosphate (NH 4 H 2 PO 4 ; 99.4% pure, Fisher Scientific) at 25 °C with a molar ratio of Ca/P = 1.67, and the amount of demineralized water is 3 times higher than that of calcite. The equation below describes the overall reaction pathway for the synthesis of Ca-HA. During the formation of Ca-HA Gel , bubbling corresponding to the release of CO 2 gas was observed. At the beginning of the reaction, the pH value increased gradually from 6 to 8 and then stabilized after 12 h of synthesis.

+ + → + + + 10CaCO 6NH H PO 2H O Ca (PO ) (OH) 6NH OH 2H O 10CO
The stabilization of the pH revealed the formation of an ammonium hydroxide buffer and the consumption of calcium carbonate. Ca-HA Powder was prepared by drying Ca-HA Gel at 105 °C for 24 h to obtain a powder with monomodal particlesize distribution. The water content calculated after drying at 105 °C was 70%. Ca-HA Gel , destined to be used for chemical characterization, was filtered and rinsed with deionized water to dissolve any unreacted phosphate. 2.3. Physical and Chemical Characterization. Powder densities were measured using a helium pycnometer (AccuPyc 1330 V2.04N). The specific surface areas of the materials (G, P, and Ca-HA) were determined on a Micromeritics Gemini 2360, by surface analysis using the Brunauer-Emmett-Teller method, and the measurements were performed by nitrogen absorption at liquid-nitrogen temperature followed by desorption (adsorption gas N 2 ; heating temperature 150 °C).

The particle-size distribution was determined by LASER granulometry using a Malvern Laser Master Sizer Hydro 2000 instrument (with ethanol as a dispersant shaken by ultrasound). The pH measurements were made after dispersion of 10 g of powder materials in 25 mL of deionized water and stirring for 1 h at room temperature, with the measurement being carried out after 2 h. The particle charge was evaluated relative to the net charge on the particles' surfaces by evaluation of the pH point zero charge (pH PZC ). This determination is based on the electrochemistry and involves the introduction of 50 cm 3 of a NaCl (0.01 M) solution in closed vials and adjustment of the initial pH to values between 2 and 12 by the addition of NaOH or HCl (0.1M). A mass of 0.15 g was added to the vials, and the final pH was measured and recorded after stirring for 48 h at room temperature. The pH PZC value corresponds to the intersection point of the titration curves obtained at one ionic strength for the suspension (pH initial ) and blank solution (pH final ).

X-ray diffraction (XRD) analysis of the powders was performed on a Siemens D5000 XRD spectrometer using Cu K α (λ = 1.540 Å) anticathode radiation (power 40 mA, 45 kV). The diffraction angle was set up from 10 to 80°2θ with a step count of 0.02°2θ and a counting time of 3 s/step. ζ-Potential measurements were performed on a Zetasizer 3000 HSA (Malvern Instruments Ltd.).

The rheological measurements were performed using a Rheostress (HAAKE RS 150 rheometer). The rheological behavior and viscosity η (Pa•s) were evaluated for shear rates ranging from 0 to 500 s -1 at constant temperature (20 ± 5 °C). The sample was placed in a cylindrical vessel, and the strain was applied by a double-helical ribbon impeller. Samples were homogenized before the shear stress τ (Pa) and viscosity versus the applied shear rate γ̇(s -1 ) were measured. The dynamic trails applied on Ca-HA Gel were performed with the same apparatus using an oscillation mode with a preselected constant shear stress of 50 Pa and a frequency ranging from 1 to 50 Hz. The loss and storage moduli measurements were recorded for 300 s at room temperature.

2.4. Formulations. Rheological tests were carried out on pure Ca-HA Gel , on binary water/gypsum (W/G), and on water/plaster (W/P) blends. Binary blends were prepared by mixing the selected material with deionized water according to different mass ratios ranging from 0.4 to 1.4. The solid concentrations of binary blends (Γ B1 ) and the W/G and W/P mass ratios were calculated as follows:

ρ ρ Γ = + = -Γ Γ = -Γ Γ V V V W G W P ; 1 ; 1 B1 s s w B G B B P B
where V s is the volume of the solid phase, V w is the volume of water, and ρ G and ρ P are the densities of gypsum and plaster, respectively. The Ca-HA Gel content considered was from 0 to 50 wt % Ca-HA Gel . The rheological behaviors of water-free formulations (Ca-HA Gel /G and Ca-HA Gel /P) with mass ratios of 50/50, 40/60, and 30/70 were not considered because the pastes obtained were very thick and almost solids). The solid concentrations of binary blends (Γ B2 ) containing either gypsum or plaster (calcium sulfate powder) with Ca-HA Gel were calculated as follows: where V s T is the total volume of the solid material equal to the sum of the Ca-HA Powder volume (V Ca-HA Powder ) and the calcium sulfate powder volume (V CS Powder ) and V w T is the total water equal to the Ca-HA water volume (V Ca-HA Water ).

Γ = + = + = - - V V V V V V V V ; ;

RESULTS AND DISCUSSION

3.1. Characteristics of the Materials. The results of the bulk density, specific surface area, and pH of the materials are given in Table 1. It can be seen that the bulk densities were used to calculate the solid concentrations in the blends. The pH of plaster is acidic, which provides an indication of the presence of phosphoric acid. The specific surface area of plaster is 3 The classification and identification of materials following the geotechnical designation reveals that plaster and Ca-HA Powder seem to have a silty-sand texture, while gypsum, having a diameter of less than 1 mm, is sandy-silt and Ca-HA Gel is identified as silty (Table 1).

XRD analysis of synthesized Ca-HA after heat treatment (calcination at 1000 °C for 10 h to crystallize the structure and eliminate calcite residue) showed the existence of a crystalline phase illustrated by hydroxyapatite in large quantities that had synthesized stoichiometrically (Ca 10 (PO 4 ) 6 OH 2 ). The XRD patterns of gypsum and plaster reveal that there are almost no common peaks. The P α and G are characterized by five major peaks at 14.7°, 25.7°, 29.6°, 31.8°, and 49.3°2θ and 11.6, 20.7, 29.0°, 31.1°, and 33.3°2θ, respectively. The presence of an anhydrite peak at 25.7°2θ appears for the first time after heat treatment (Figure 1). The difference in the structures of both plaster and gypsum is related to the heat treatment, which modifies the crystalline lattice. The study carried out by Mandal et al. demonstrated that heating the gypsum (CaSO 4 •2H 2 O; 98% pure) at 90 and 350 °C for 10 h revealed the presence of two kinds of water: loosely and strongly held water molecules. [START_REF] Mandal | Anion water in gypsum (CaSO 4 •2H 2 O) and hemihydrates (CaSO 4 • 1 / 2 H 2 O)[END_REF] Consequently, the loss of water molecules modifies the monoclinic structural morphology of gypsum, where the elimination of 1.5 water molecules allows the apparition of hexagonal hemihydrates. Above 200 °C, the symmetry changes to the orthogonal configuration due to the elimination of strongly held water molecules. As shown in Figure 2, the shear stress increases with the shear rate when the stress exceeds a threshold value τ 0 (1.425 Pa), indicating that Ca-HA Gel exhibits a shear thinning or, specifically, a viscoplastic behavior. Studies published by Knowles et al. carried out on a commercial Ca-HA suspension agree with this obtained result. [START_REF] Knowles | Characterisation of the rheological properties and zeta potential of a range of hydroxyapatite powders[END_REF] The viscosity reaches an apparent value of η ∞ = 0.144 Pa•s above γ̇= 200 s -1 . The critical shear rate corresponding to the Ca-HA dispersion is characterized by an extrapolated maximum yield stress τ 0e (2.289 Pa) and a minimum yield stress τ 0 (1.425 Pa), related to the response of the gel to the deformation. The viscosity of Ca-HA Gel is 0.144 Pa•s at a higher shear stress, and above this value, the suspension flows. The increase in the shear stress leads to the disruption of the gel structure, and when the shear stress stops, the suspension tends to recover its original structure (i.e., the flow is reversible).

RHEOLOGICAL BEHAVIOR

The viscosity does not only depend on the shear rate, but it is also correlated to the duration of the applied stress. The internal structure of Ca-HA Gel undergoes an initial homogenization at low shear rate values corresponding to low levels of plasticity. The constraint force is too small to disorganize the structure; it requires an increase in the magnitude of the applied force for the flow to occur. Beyond the shear stress threshold, with an increasing shear rate, the flow is enhanced by the gradual dispersion of particles, promoting the separation effect, which causes a decrease in the viscosity. The viscosity (apparent viscosity) value stabilizes at high shear rates. Once the apparent viscosity is reached, there is no further modification observed for the structure of Ca-HA Gel .

Dynamic Rheological Properties of Ca-HA.

To explain the viscoelastic behavior observed for Ca-HA Gel , dynamic oscillatory tests were carried out to investigate the internal structure and response of the gel during Ca-HA synthesis. The dynamic rheological (oscillatory) test aims to determine the storage (G′) and loss (G″) moduli. Figure 3 represents the results of the frequency-sweep response of a Ca-HA suspension after 12, 24, 48, and 72 h of synthesis.

During synthesis, the storage and loss moduli (G′ and G″) rose with an increase in the frequency and their values decreased with the aging of the gel.

At low frequencies, the 12 h age dispersion (Figure 3a) seems to have a heterogeneous structure at the beginning of the periodic oscillation application. It is observed that G′ is greater than G″, which corresponds to an elastic behavior. Therefore, when the frequency is increased, the modulus values increase strongly because of the dispersed state of the particles.

For a 24 h age dispersion (Figure 3b), the values of the storage and loss moduli are almost zero below 15 Hz. The strain response at constant stress after 24 h of synthesis is delayed compared to the 12 h age dispersion measurement. This change in behavior, identified at low frequencies, is due to the evolution of the reaction. This response is related to precipitation of the Ca-HA suspension occurring when the pH value buffers around 8, leading to the formation of a colloidal suspension. The behavior remains elastic at high-frequency values.

After 48 h of synthesis (Figure 3c), we observe a complete change in behavior at low frequencies. For frequencies below 25 Hz, G′ and G″ increase exponentially at low stress. Below 25 Hz, G″ is greater than G′, indicating viscous behavior. This simply means that the suspension follows a constant structural deformation where the stress is proportional to the strain rate. When the frequency reaches values higher than 25 Hz, G″ becomes smaller than G′, indicating an elastic behavior: the stress is proportional to the deformation. However, the viscosity increases progressively, indicating the agglomeration of colloids and the formation of gel-paste calcium phosphate. This result is important because it shows the transition between a suspension and pasty behavior. As demonstrated by Knowles et al. [START_REF] Knowles | Characterisation of the rheological properties and zeta potential of a range of hydroxyapatite powders[END_REF] and Bao et al., [START_REF] Bao | Rheological behavior of aqueous suspensions of hydroxyapatite (HAP)[END_REF] the mutual approach of particles in solution is correlated to the positive charge of the Ca-HA surface. The behavioral changes show the formation of gel particles or gel-paste structures after 2 days of reaction, exhibiting a viscoelastic behavior. The dynamic rheological property study carried out by Liu et al. [START_REF] Liu | Rheological properties of concentrated aqueous injectable calcium phosphate cement slurry[END_REF] on a calcium phosphate cement slurry (a composite of equimolar tetracalcium and dicalcium phosphate anhydrous) under the conditions of a 2.0 powder/liquid ratio and a 1% strain showed that, at the range of low strain and low frequency, the slurry changed its behavior with frequency increases. The slurry began to flow and presented a shear-thinning behavior.

For the 72 h age dispersion (Figure 3d), the loss modulus values are higher than the storage modulus values, thus revealing a viscous behavior. Nevertheless, when the frequency increases and the strain exceeds the yield strain, the viscosity decreases. The decrease in the G′ and G″ values is due to gel aging after 48 h of synthesis, which causes deformation of the gel structure. The disruption of extractible water molecules caused by bonds between particles leads to flow and demonstrates that Ca-HA Gel exhibits a shear-thinning behavior.

4.3. Rheological Behavior of Binary Blends. The viscosity of the binary matrix was measured and plotted versus the concentration of added powder, as illustrated in Figures 4 and5. W/G and W/P blends exhibit shear-thinning behavior. Figure 4 shows that the viscosity of W/G blends decreases with decreasing gypsum concentration until an almost constant apparent viscosity for W/G ratios above 0.6 was achieved. The viscosity plateau reveals no influence of water on the hydrated gypsum particles. Hydrated gypsum particles present similar double-electric layers of a common charge (negatively The description of the behavior of Ca-HA Gel is based on two conditions: if G′ > G″, the behavior is elastic; if G′ < G″, the behavior is viscous. charged), which generates electrostatic repulsion between particle surfaces (due to the overlap of the electrical double-layer theory). Gypsum particles are mainly positively charged and are covered by a layer of positive average charge. This positivecharge layer stabilizes the dispersion through particle repulsion, avoiding particle aggregation. At rest, particles tend to precipitate. If a stress is applied, the suspension is dispersed, leading to a decrease in the viscosity. The gypsum concentration controls the dispersed state; the repulsive forces tend to stabilize the suspension because the particles are negatively charged. The apparent viscosity is obtained at a high shear rate, where gravity forces are negligible. It can be seen that, above a 0.6 W/G ratio, the interparticle distance is large enough to make the apparent viscosity independent of the gypsum content.

With regards to the plaster blends (W/P), the viscosity decreases with decreasing plaster concentration. The hydration of hemihydrate particles, requiring 1.5 water molecules, leads to crystallization. In general, the contact of plaster particles with water triggers dissolution, which induces the saturation or supersaturation of the solution in Ca 2+ and SO 4 2-, leading to the formation of gypsum nuclei. [START_REF] Singh | Calcium sulphate hemihydrate hydration leading to gypsum crystallization[END_REF] Consequently, the hydration control scenario is correlated with the nucleation and crystalgrowth mechanism. [START_REF] Bosbach | Gypsum growth in the presence of background electrolytes studied by Scanning Force Microscopy[END_REF][START_REF] Sievert | Hydration of anhydrite of gypsum (CaSO 4 .II) in a ball mill[END_REF] The high specific surface area of plaster particles plays a major role in the multiplication of intergranular links. Interaction between plaster particles occurs by electrostatic forces (i.e., the spread of ions in an extended diffuse layer) and attractive forces (van der Waals) responsible for adhesion. Furthermore, adhesion is described as the formation of two Ca 2+ /SO 4 2-layers bound together by ionic interaction, constructing a growth slice. [START_REF] Singh | Calcium sulphate hemihydrate hydration leading to gypsum crystallization[END_REF] The predominance of reactive sites on the grains attracts cations, which form a rigid network that surrounds the layer of adsorbed water molecules. It is also known that the attraction of plaster particles requires cohesion between atoms to form the contiguous flakes. The development of this structure (agglomerated structure or plaster paste) depends on the amount of hydration water. According to the stochiometry, the α-hemihydrates require a W/P ratio higher than 0.186 to be transformed into gypsum. To improve the consistency of the blends, the hydration rate was chosen from a range of 0.4-1.4. In the dynamical study in the linear viscoelastic domain published by Finot et al., the time-sweep test at 0.015% (ω = 10 rad•s -1 ), carried out on W/P = 0.7, indicates that the rise in the storage modulus reflects the hydration and setting processes. This reveals an increase in the level of crystallized gypsum. [START_REF] Finot | Correlation surface forces with surface reactivity of gypsum crystals by atomic force microscopy. Comparison with rheological properties of plaster[END_REF] Obviously, the viscosity decreases steadily with an increasing amount of water, and the application of shear stress destroys the links between particles, facilitating decoagulation of the structure and leading to flow.

Plaster hydrates during viscosity measurement under stress. The plaster viscosity value depends on the volume fraction of the solids (solid concentration Γ B1 ), particle size, and water excess (free of water after full plaster hydration). At a high shear rate, it is possible to measure the apparent viscosity of the dispersion of agglomerated hydrated particles. The apparent viscosities of W/G blends are lower than those of W/P blends. Because plaster hydration held the measurement long, the difference observed in the apparent viscosity between W/G and W/P cannot be solely attributed to hydration. The difference observed in a specific surface area may provide an explanation. Figure 5 shows that the viscosity increases when the solid concentration (Γ B2 ) increases, with the latter being calculated relative to the water (by wt %) present in Ca-HA Gel . It should be noted that the formulations having Ca-HA Gel /G or Ca-HA Gel /P ratios higher than 1.5 are solids and their rheological properties were not characterized. The water content in Ca-HA is more than sufficient to hydrate calcium sulfate hemihydrate particles. In this case, an almost solid paste is formed. Replacing 20% of Ca-HA Gel by gypsum or plaster produced little influence on the viscosity, indicating that for this blend there is little influence of calcium sulfate on the rheological properties of Ca-HA Gel . With regards to the already hydrated gypsum particles, increasing the gypsum concentration of the blend up to 40% produced a slight increase in the viscosity. Blends containing more than 40% of gypsum are almost solids. The addition of plaster to Ca-HA requires a high stress to cause the flow. Plaster particles in contact with Ca-HA water cause plaster hydration, creating strong bonds that encompass the Ca-HA particles and make the flow of the blend difficult. Hydrated plaster strongly increases the viscosity of Ca-HA Gel by its own weight. Figure 6 shows the thixotropic behavior of binary blends containing a mass ratio blend of 80% of Ca-HA and 20% of plaster and a mass ratio of 80% of Ca-HA and 20% of gypsum, indicating that the blends recover their initial structure. In this case, viscosity values increase slightly compared to the viscosity of Ca-HA Gel while keeping their flow characteristics.

4.4. Interaction between Particles. The interpretation of interactions between these compounds seems to be complex. The ζ potential was considered as a means to understanding the interactions, based on the net charge of the particles, and hence describing the suspension structure. The ζ potential indicates the repulsion forces between particles promoting the dispersion stability. Figure 7 shows that the Ca-HA particles are characterized by a negative ζ potential, indicating that the particles are negatively charged (≈-20 mV) [START_REF] Knowles | Characterisation of the rheological properties and zeta potential of a range of hydroxyapatite powders[END_REF] and the gypsum particles are positively charged (≈20 mV). [START_REF] Peng | Adsorption characteristics of water-reducton agents on gypsum surface and its effects on the rheology of gypsum plaster[END_REF] The Ca-HA ζpotential value indicates the instability of the gelled particles. Adding gypsum to Ca-HA leads to charge modification. The presence of repulsive forces is correlated to the high value of the ζ potential predicting the stability. These results support the assumption of sulfate attraction on the Ca-HA surface, indicating dispersion of the particles. However, sulfate adsorption changed the charge distribution of Ca-HA, generating the change from negative to positive charge. The adsorbed amount of anionic sulfate on the Ca-HA particles is equal to the Ca-HA electrical charge sites, and the coating of gypsum sulfate Ca-HA particles is the main reason for modification of the ζ potential. The surfaces of the plaster particles are characterized by a positive charge, and the formation of a plaster solid body is explained by cohesion between particles and illustrated by a low ζ-potential value.

Figure 7 represents the ζ potential of formulations containing Ca-HA and different proportions of the W/P ratio. With an increase in the amount of plaster, the ζ potential remains constant (with values similar to those of the plaster at the beginning of the hydration process) and the surface charge of Ca-HA changes. The nucleation process of hemihydrated particles starts by the formation of Ca 2+ /SO 4 2--bound double layers, and the concentration of sulfate decreases, indicating that it is consumed during hydration. However, fixation of the remaining sulfate in solution on the Ca-HA surface leads to a change in the charge.

CONCLUSION

The rheological test carried out reveals shear thinning and the thixotropic behavior of Ca-HA Gel for formulations considered. The internal structure of Ca-HA Gel exhibits viscoelastic behavior after 48 h of synthesis.

The positively charged gypsum particles cause repulsion of the particles (repulsive forces), and the viscosity of the W/G (0.2-1.4) blends remains constant. Plaster is characterized by attractive forces, and the particles are positively charged. The hydration of plaster leads to cohesion between particles, which explains the increase in the viscosity of the W/P blends compared to the W/G blends. The ζ potential of Ca-HA is negative, and dissolution of the calcium sulfate particles saturates the solution in sulfate ions (SO 4 2-); the charge of the particles becomes positive because of fixation of SO 4 2-on the Ca-HA surface. The adsorption of sulfate on the Ca-HA particle surfaces influences the structures of the formulations and explains why the viscosity is independent of the Ca-HA Gel content. Plaster hydration plays a major role in the Ca-HA Gel stabilization structure, and gypsum particles undoubtedly enhance the hydraulic performance. Modification in the ζ potential can affect the retention capacity, reactivity, and mechanical resistance. These aspects will be the subject of a later study.

  times higher than that of gypsum. The specific surface area of Ca-HA Powder indicates the surface reactivity of the particles. The net charge of the materials can be evaluated by referring to the pH of the aqueous solutions with the metal surface in a neutral electric state and pH PZC . The plaster and Ca-HA particles are positively charged (pH < pH PZC ), and the gypsum particles are negatively charged (pH > pH PZC ).

4. 1 .

 1 Steady Rheological Properties of Ca-HA Gel . Apatite gel is characterized by a specific texture caused by absorption of water and interaction between fine particles. The identification of the behavior and the structure Ca-HA Gel requires an understanding of the influence of parameters such as material deformation, strength, and stability.

Figure 1 .

 1 Figure 1. Powder XRD pattern of Ca-HA heated for 10 h at 1000 °C [hydroxylapatite (+), hydroxyapatite (×), plaster (∧), gypsum (*), and anhydrite (○)].

Figure 2 .

 2 Figure 2. Rheological behavior and viscosity of Ca-HA Gel after 48 h of synthesis: (a) τ = f(γ); (b) log(τ) = f(γ) [-×-, shear stress (Pa); -+ -, viscosity (Pa•s)].

Figure 3 .

 3 Figure 3. Frequency-sweep curves of Ca-HA as a function of the synthesis time (-Δ-, storage modulus G′, -□ -, loss modulus G″).The description of the behavior of Ca-HA Gel is based on two conditions: if G′ > G″, the behavior is elastic; if G′ < G″, the behavior is viscous.

Figure 4 .

 4 Figure 4. Viscosities of W/G and W/P blends at different mass ratios.

Figure 5 .

 5 Figure 5. Viscosities of Ca-HA/G and Ca-HA/P formulations.

Figure 6 .

 6 Figure 6. Example of a thixotropic hysteresis loop of composites 80% Ca-HA/20% G and 80% Ca-HA/20% P [-*-, Ca-HA Gel (η ∞ = 0.144 Pa•s); -×-, 80% Ca-HA Gel /20% P (Γ = 0.369; η ∞ = 0.348 Pa•s); -+-, 80% Ca-HA Gel /20% G (Γ = 0.384; η ∞ = 0.318 Pa•s)].

Figure 7 .

 7 Figure 7. ζ potential of Ca-HA/G and Ca-HA/P formulations without added water [--, ζ potential of gypsum (≈20 mV); -• -, ζ potential of plaster (≈5 mV); -•• -, ζ potential of Ca-HA (≈-20 mV); -□ -, ζ potential of Ca-HA/G formulations; -Δ-, ζ potential of Ca-HA/P formulations].

  2.1. Calcium Sulfates. The gypsum (G, CaSO 4 •2H 2 O) and the plaster calcium sulfate (P α , CaSO 4 •0.5H

  2

Table 1 .

 1 Physical Characteristics of the Materials

	G Powder P Powder	Ca-HA Powder	Ca-HA Gel

This study has led to the understanding of the rheological and physicochemical behavior on water/calcium sulfates and hydroxyapatite/calcium sulfate blends. These results will be useful in the formulation of a permeable reactive barrier.

■ AUTHOR INFORMATION

Corresponding Author *Tel: +(00)33563493241. E-mail: Mohamed.Raii@mines-albi. fr.

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

The authors are thankful to PRAYON Belgium for providing gypsum and plaster and their interest in the present study.


Julien Lardy

Thomas Robert

Xuguang Wang

Proposition of an algorithm

come

Proposition of an algorithm for simulating natural variability in arm grasping posture when rotating a spherical object Introduction Digital human models are being more and more used in early phase of a product/workplace design. However, human motion simulation still faces the well-known problem of motor redundancy, implying that humans have more degrees of freedom than necessary to perform a given task. As a consequence, humans have an infinite number of possibilities to do one task. Solutions used nowadays in motion simulation are mainly based on optimization trying to find a unique solution that satisfies a set of constraints while minimizing an objective function (see [1] for review). Nevertheless, some authors suggest that redundancy should be considered as an advantage rather than a problem to be solved. Indeed, Latash considers the human redundancy as an ability to adapt to the environment [START_REF] Latash | Fundamentals of Motor Control[END_REF]. Moreover, studies showed that different ways to perform a task can coexist, questioning the pertinence of looking for a unique solution. Then, natural motion variability should also be simulated. Grasping and manipulating an object with the hand is one of the routine tasks of human beings and is frequently encountered in workplace evaluation and product design. In case of rotating a spherical object, a recent study showed a large variation in arm posture when grasping the object for a given rotation amplitude [START_REF] Lardy | [END_REF]. The work presented here aims at proposing a stochastic algorithm to predict the grasping posture as a function of the intended manipulation.

Methods

The rotation task consists in grasping and rotating a sphere. Rotation amplitudes range from -315° in counterclockwise direction (CCW) to 270° in clockwise direction (CW). Actual subjects performed this task in order to compare our simulation results to real data.. Experimental results and details of the protocol have been published in [START_REF] Lardy | [END_REF]. Hand palm is considered as a rigid plan formed by the head of the second metacarpal bone and the wrist styloids. Although there are multiple contact points between the palm and the sphere, we choose to limit our analysis to a unique contact point on the sphere, where the sphere and the palm share the same normal vector. From the shoulder to the palm, seven degrees of freedom are defined (3 for the shoulder, 2 for the elbow and 2 for the wrist). Joint limits avoidance is based on Wang's work [4] which takes into account coupling effect of joint axes allowing us to have more natural simulated postures. Shoulder and wrist circumduction cones are based on data from the literature ( [4] for the wrist and [6] for the shoulder), while axial limits of the upper arm are from [5]. The definition of the grasping posture is done as follow: 1-estimation of all the possible grasping points allowing an amplitude of rotation required by task; 2-a uniform probability distribution is attributed to these admissible points; 3-a grasping point is selected randomly; 4-the palm orientation around the selected grasping point is defined in order to maximize the shoulder-wrist pivoting angle margin at the beginning of the rotation, given the constraint that it must be inferior or equal to the pivoting margin at the end of the rotation; 4bis-if this constraint cannot be satisfied, another point is selected (step 3); 5-the elbow position is compute in the middle of the pivoting margin once the palm orientation is determined.

Results

Simulations were based on one subject's anthropometric data. Due to simulation constraints (rigid hand and no shoulder motion), we were only able to simulate grasping postures for the rotation amplitudes from -315° to 270°. Twenty simulations were performed for each amplitude. Figure represents the forearm (prono-supination) and shoulder-wrist (pivoting) angles for simulated and real postures. Globally, one can observe a similar evolution between the simulated and actual angles. Also, we have to note a tendency to have higher values for simulated pivoting angle in counter-clockwise rotations, i.e. a lower elbow position. Finally, as for the real data, simulated pivoting margins at the beginning and end of rotation show same behavior described in one of our previous work [START_REF] Lardy | [END_REF].

LARDY, Julien, ROBERT, Thomas, WANG, Xuguang

Conclusion

First conclusion of this work is the ability of our algorithm to simulate arm grasping postures in agreement with experimental data. Moreover, despite a quite simple probability distribution, simulations show interesting results in term motion of variability. Nevertheless this model could be improved taking into account individual joint limits. Finally, other probability distributions should also be tested.

Biblio